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Abstract. Diffusion preserves the positivity of concentrations, therefore, multicomponent diffu- 
sion should be nonlinear if there exist non-diagonal terms. The vast variety of nonlinear multicom- 
ponent diffusion equations should be ordered and special tools are needed to provide the systematic 
construction of the nonlinear diffusion equations for multicomponent mixtures with significant in- 
teraction between components. We develop an approach to nonlinear multicomponent diffusion 
based on the idea of the reaction mechanism borrowed from chemical kinetics. 

Chemical kinetics gave rise to very seminal tools for the modeling of processes. This is the 
stoichiometric algebra supplemented by the simple kinetic law. The results of this invention are 
now applied in many areas of science, from particle physics to sociology. In our work we extend 
the area of applications onto nonlinear multicomponent diffusion. 

We demonstrate, how the mechanism based approach to multicomponent diffusion can be in- 
cluded into the general thermodynamic framework, and prove the corresponding dissipation in- 
equalities. To satisfy thermodynamic restrictions, the kinetic law of an elementary process cannot 
have an arbitrary form. For the general kinetic law (the generalized Mass Action Law), additional 
conditions are proved. The cell-jump formalism gives an intuitively clear representation of the 
elementary transport processes and, at the same time, produces kinetic finite elements, a tool for 
numerical simulation. 
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1. Introduction 

1.1. Linear Diffusion: from Graham and Fick 
to Einstein, Onsager and Teorell 

1.1.1. Fick's Law 

The first prominent equation of diffusion is Fick's Law. According to this law, the diffusion flux J 
is proportional to the antigradient of the concentration c: 

J=-Dgradc. (1.1) 

The time derivative of the concentration is the negative of the divergence of the flux: 

^ = -dwJ = DAc, (1.2) 

where A is the Laplace operator. 

This statement is closely related to the Gauss-Ostrogradskii theorem 

(divJ) dV = jjj-ndS. (1.3) 

v s 

The left side is an integral over the volume V, the right side is the surface integral over the boundary 
S of the volume V, S = dV, and n is the outward pointing unit normal field of S. The right-hand 
side represents the total flow across the boundary "out of the volume V". The theorem was first 
discovered by J. L. Lagrange in 1762, then later independently rediscovered by C. F. Gauss in 
1813, by G. Green in 1825 and in 1831 by M. V. Ostrogradsky. According to this theorem, the 
diffusion equation d t c = — div J is just a conservation law: all changes of concentration are caused 
by the flux only. 

In his work on diffusion law, A. Fick used the conservation of matter and the analogy between 
diffusion and Fouriers law for heat conduction (1822), or Ohms law for electricity (1827). De- 
velopment of the fundamental law of diffusion was inspired by the Graham's investigations on the 
diffusion of salts in water 115511 . in which he studied and compared the diffusibility of different salts. 

Before his study of diffusion in liquids, Graham studied diffusion in gases (1833). In 1863, J. 
C. Maxwell calculated the diffusion coefficients in gases from the Graham data. The results are 
amazing: "His coefficient of diffusion of C0 2 in air is accurate at ±5%. Isn't it extraordinary?" 

Maxwell's theory of diffusion was based on gas kinetics and mean fee path estimates. 




1.1.2. Einstein's Mobility 

In his theory of Brownian motion, A. Einstein ll24l developed the microscopic theory of the diffu- 
sion coefficient for diluted particles in a liquid and compared two processes: the motion of particles 
in a liquid under a constant external force K , and diffusion. For a given K, each particle has the 
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average velocity mK where the coefficient m characterizes mobility of particles. (We use m for 
mobility and reserve [i for chemical potential.) For spherical particles in liquid, 

D7rr]r 

where rj is the viscosity of the liquid, r is the radius of particles, and 67rr]r is the Stokes friction 
force. 

This approach results in a very useful relation for the diffusion coefficient: 

RT 

D = m— = mk B T, (1.5) 

l\ A 

where R is the gas constant, Na is the Avogadro constant, k B is the Boltzmann constant. The 
coefficient m is called mobility or the Einstein mobility. 

Graham's experimental research was extended to solids by W. C. Roberts-Austen [91]. He used 
Fick's equation to determine the diffusion coefficient ll79l . In 1922, S. Dushman and I. Langmuir 
j[22| proposed to use the Arrhenius law for diffusion coefficient: 

D = D exp {—QJ kT) , (1.6) 

where Q is a constant, which we now recognize as the activation Gibbs energy of diffusion AG. 
More precisely, AG includes two terms: AG = AH — TAS, where AH is the activation enthalpy 
and A S is the activation entropy. 

They checked this law by their own experiments with the diffusion of thorium through tungsten 
and found satisfactory agreement. Even better agreement was found with the published results of 
W. C. Roberts-Austen's experiments. 



1.1.3. Teorell Formula 

The mobility-based approach was further applied by T. Teorell H101II . In 1935, he studied the dif- 
fusion of ions through a membrane (Fig.Q]). He considered a system of an ideally dilute solution of 
binary univalent strong electrolysis at the same temperature in water. The boundary is considered 
to be a membrane with strong a electrolyte in the presence of water. The solutions are assumed to 
be kept homogeneous on both sides of the membrane up to the boundary by some form of convec- 
tion. The ionic mobilities within the membrane are assumed constant and may be different, and 
the membrane is not permeable for water. Heat effects, special membrane effects and chemical 
reactions are ignored. 

He formulated the essence of his approach in the formula: 
the flux is equal to mobility x concentration x force per gram ion. 

This is the so-called Teorell formula. 

The force consists of two parts: 

1. Diffusion force caused by concentration gradient: —RT-j^. 
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Figure 1: Teorell's model HIOIH . A system of ideally dilute solutions in water of binary univalent 
strong electrolytes DA, M'B', M"B" is considered. The boundary membrane is permeable 
for all ions (cations D, M,... and anions A, £>,...) but the movement of water is prevented. The 
ionic concentrations outside are kept constant. Also Df is maintained constant. Accordingly, DA 
is a steadily diffusing electrolyte. No other electric field is present besides that due to diffusion 
potential. No chemical reaction takes place. The steady state of a system of this nature with a 
steady diffusion is characterized by a constant ratio series: Tr = Tr* = ''' = ~w = Tr* = '''- 

2. Electrostatic force caused by electric potential gradient: q^. 

Here R is the gas constant, T is the absolute temperature, c is the concentration, q is the charge 
and if is the electric potential. 

In these notations, the Teorell formula for the flux J is 

/ RTdc d<p\ 

J = mc[ -r + ^TT 1 ( L7 ) 

\ c dx dx J 

(m denotes mobility; here we slightly modernize notations). It may be worthwhile to introduce the 
reference equilibrium concentrations vector c* and write the diffusion force in the form 

RTdc dln(c/c*) 
— = -RT • (1-8) 

c dx dx 

This expression allowed Teorell to find the concentration jump and the electric potential across the 
membrane caused by the joint action of diffusion and the electric field, when mobilities of various 
components are different. 

1.1.4. Onsager's Linear Phenomenology 

In 1931, L. Onsager ll83l |84| included diffusion in the general context of linear non-equilibrium 
thermodynamics. For multi-component diffusion, 



J i = Y, L n x ji (1-9) 
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where Jj is the flux of the ith component and Xj is the jth thermodynamic force (for pure dif- 
fusion, this is the space antigradient of the jth chemical potential divided by T). After lineariza- 
tion near equilibrium, this approach gives for perfect systems (for which the chemical potential is 

RT\n(c/c*)): 

1 , R 



Xj = -— gradfij- = — -grades 
1 c j 

Ji = -^2^— grade,-; 

3 j 

= - divJ * = R ^2 L a-^ 

i 3 



(1.10) 



J 



where c* are equilibrium constants (c* is the point of linearization), deviations of Cj from c* are 
assumed to be small, A is the Laplace operator, and Lij = Lji is the matrix of the coefficients. Its 
symmetry follows from microreversibility. 

The system (11.10b has one attractive property. Let us consider this system in a bounded domain 
V with smooth boundary and with zero fluxes through its boundary: (n, grade,) = at any point 
of dV at any time (n is the vector of the outer normal). The positive quadratic functional 

Z i Jv Ci 

is the second-order approximation to the relative entropy (or the so-called Kulback-Leubler diver- 
gence, see the review paper H71 ) 

Skl = Y,L< x )^( 2 ^t) dx - (L12) 



V 



Let us calculate the time derivative of S2 due to the system (11.10b . Using the Gauss-Ostrogradskii 
formula dl.31) we get for the positive semidefinite matrix L: 

dS 2 f v-^ fVci VcA 

ir = -U£MT'^) dI - ' (U3) 

where ^ j is the inner product of the space vectors. 

Therefore, ^ < if the symmetric coefficient matrix is positive semidefinite (this means 
that for any vector £ the following inequality holds: £\ . L^^j > 0). 

The Onsager form of the diffusion equations is correct near the equilibrium but violates the 
obvious physical requirement: the diffusion flux of the ith component is zero if its concentration 
has zero value: the flux vanishes with the concentration. The Teorell formula satisfies this require- 
ment. Fick's law also satisfies this requirement in the following sense: if for nonnegative smooth 
c(x) the concentration vanishes at some points, then at these points the flux vanishes too (because 
these points are minimizers of concentration and the gradient vanishes there). 
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For isotropic non-perfect systems we have to use the generalized thermodynamic forces in 
Onsager's form of the diffusion law: 



v - d f 



grade.,-, (1-14) 



where T/(c, T) is the free energy density. 

Let us denote $^ = (d 2 f / dcidcj) c=c *. In this notation, 



dt 

The quadratic form 



Xj = - ^ $ ife gradc fc ; J { = ^ LyXj = - ^ I ^ ^$ jfc J gradc fc ; 

k j k \ j / 

-divj, = (j2 L ij^ At 



(1.15) 

dci 

Ck ■ 



F 2 = ljYl ^{Cj ~ C*)( Cfc - C*) dx 



is positive definite because F is strictly convex. For positive definite L, F 2 decreases in time due 
to diffusion. Indeed, in a bounded domain V with a smooth boundary and without fluxes through 
the boundary we get analogously to (11.131) : 



dF 2 
dt 



/X (Z^ Vcfc j L « (E^ Vc ^ dx ^ °- (L16) 



For non-isotropic diffusion (for example, in crystals), the coefficients L have two pairs of indexes: 
Liaj/3, where i,j correspond to components and a, (3 correspond to the space coordinates. The 
forces and fluxes also have these two indexes and 

Jia — L ia jpXjp . 

In all cases, the diffusion equations in Onsager's form do not describe the non-diagonal terms 
(the influence of gradients q on fluxes of Cj for i ^ j) properly near zeros of concentrations. These 
equations are applicable near a reference point c* > only. 

Non-diagonal diffusion must be non-linear. This simple remark is so important that we will 
explain it in detail. Let diffusion be non-diagonal and linear: 



d t Ci = ^ DijAcj . 



Assume that D\ 2 ^ and consider the state with c 2 — . . . — c n — 0. At this state, 

d t c 2 = D u Aci . 

If Di 2 Aci(x) < at some points then c 2 (x) becomes negative at these point in a short time. 
Therefore, linear non-diagonal diffusion does not preserve positivity of concentrations. 
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1.2. Mechanisms of Nonlinear Diffusion 
1.2.1. Jumps on the Surface 

In 1980, A.N, Gorban, V.I. Bykov and GS. Yablonskii H51 proposed a model for diffusion in 
monolayers of reagents on the surface of a catalyst, which is based on the jumps of the reagents on 
the nearest free places. This model was used for CO on Pt oxidation under low gas pressure. 

The system includes several reagents Ai, A 2 , ■ ■ ■ A n on the surface. Their surface concentra- 
tions are ci, c 2 , . . . c n . The surface is a lattice of the adsorbtion places. Each reagent molecule 
fills a place on the surface. Some of the places are free. We use Z = A for a free place and the 
concentration is z = cq. The sum of all c, (including free places) is constant: 



i=0 



const . 



The jump model gives for the diffusion flux of A* (i = 1, . . . , 

Ji = -Di[zVci - qVz] . 
Therefore, the corresponding diffusion equation is: 

-div Ji = Di[zAci — CiAz 
Due to the conservation law, 

n 

z = b - ^ °i » 
i=l 

and we have the system of n diffusion equations: 



n): 



dcj 
dt 



(1.17) 



(1.18) 



Ji = 

dcj 
dt 



-Di 



Di 



i=l J \i=l 
n \ / n 

b - Oj j Acj + CjA I Ci 



1=1 



J=l 



(1.19) 



It is straightforward to check that when c, > for all x then d t Ci > for Cj = 0. This is a necessary 
condition for preservation of positivity. 

If we assume that all particles can exchange their positions with their closest neighbors then a 
simple generalization of (11.171) . ( 11.191) appears: 



Ji = ~ DijlcjVci - CiVcj] ; 



dc 



(1.20) 
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where Dy = Dji > is a symmetric matrix of coefficients which characterize the intensities of 
jumps. 

The entropic Lyapunov functional for (11.191 ), (11.201 ) has a simple traditional form of perfect 
relative entropy: for any reference vector of concentrations c* (c* > 0) 

Skl = j E c * ln (jl) dx - (L21) 

Remark: the free place entropy should be obligatorily included into Skl- 

Simple algebra gives that in a bounded domain V with a smooth boundary and without fluxes 
through boundaries 

1f-p4te Vc '-H 2dl - - a22) 

This inequality provides the Lyapunov stability of diffusion. 

It is worth mentioning that the thermodynamic inequality (11.221) requires only the non-negativity 
of coefficients Dij and does not imply any requirements on the matrix fasa whole (like positive 
definiteness). Another form of the thermodynamic inequality makes the formula for the entropy 
production more transparent: 

dS KL _ ^ n f f„,fCj\ „, /C,; X 



dt 



Y, A; f y (qV In (|) - Cj V In (^J) dx < . (1.23) 



This system of models was further developed by A.N. Gorban and H.P. Sargsyan [53] and 
published in a book 11461 . 



1.2.2. Diffusion in Solids as Reaction: from Frenkel to Eyring 

The physical idea of the quasi-chemical representation of diffusion in solids belongs to Yakov 
Frenkel Il36ll37ll . He introduced both the vacancy and the interstitial mechanisms of diffusion and 
found some rate constants from experimental data. 

Thirty years later, F. C. Frank and D. Turnbull developed the Frenkel theory further [|34ll . They 
studied the diffusion of copper in germanium. This diffusivity is very rapid. They proposed that 
the copper could be dissolved in two states, interstitial and substitutional. For the interstitial state 
the solubility of copper is two orders of magnitude less and the diffusivity many orders of mag- 
nitude greater than in the substitutional state. The conversion of these states is effected by lattice 
vacancies. 

The quasi-chemical theory of diffusion and viscosity was developed also by H. Eyring with 
co-authors ll66l . Eyring developed the theory of absolute reaction rates for chemical reactions in 
gases ll28l and in condensed phase H116H and then applied these ideas to transport phenomena. 

In this theory, the transport process is represented by an ensemble of elementary events. Each 
elementary event is represented by the creation or disintegration of an activated complex. The rate 
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of the elementary process is given by the concentration of activated complexes, multiplied by the 
rate at which they decompose. 

The main constructive hypothesis is that it is possible to calculate the concentration of acti- 
vated complexes by equilibrium statistical thermodynamics: the complex concentration is in quasi- 
equilibrium with the stable components. Each complex has its "internal translational" degree of 
freedom. On the surface of potential energy this corresponds to the "reaction path". Complexes 
move along this path. The velocity of this motion is assumed to be just a thermal velocity and is 
proportional to y/T. 

The additional reaction path degree of freedom has its own kinetic energy and, therefore, in- 
creases the complex heat capacity. We have to take this into account in the calculation of the 
equilibrium constant. 

Collective models of diffusion were proposed too. One of the earliest collective model is the 
Z. Jeffries "ring mechanism" with 4 or more atoms. More on the history of solid-state diffusion is 
presented in the review HI 0411 and in a modern textbook 11781 . 

On the surface, there are various mechanisms for collective diffusion ll88l as well. Elementary 
events for these mechanisms involve many atoms simultaneously. A dynamic description of non- 
linear multicomponent diffusion requires a unified framework that should satisfy basic physical 
principles. 

1.2.3. Ginzburg-Landau Free energy and Cahn-Hilliard equation 

The processes of phase separation has remained for a long time an important source of problems 
and ideas for the theory of nonlinear diffusion. The analogue for Fick's equation is the Cahn- 
Hilliard equation |[T2l . 

The Cahn-Hilliard equation in its simplest form has the standard Onsager form, the flux is 
proportional to the force, the force is the gradient of the chemical potential: 



If we compare this equation to the Teorell formula then we immediately find the missed factor 
c (concentration). We will return to the problem of the proper prefactor in the Cahn-Hilliard 
equations later. The main specificity of the Cahn-Hilliard equations is the form of the free energy 
and the chemical potential lfT2~l[TTTl (the Ginzburg-Landau form): 



The term 7(Vc) 2 in the free energy penalizes ovr sharp gradients and, in particular, models the 
interface energy. 

According to (11.241) and (11.251) . the Cahn-Hilliard equation reads 



J = —DV/jL , — = -divJ = DAfi . 



(1.24) 




(1.25) 




(1.26) 
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If the chemical part of the free energy, f c (c) is not convex then phase separation is possible. If at 
the point c this function is concave (spinodal) then without the term with A 2 the constant solution 
c(x) = c for the diffusion equation becomes unstable to any perturbation (negative diffusion coef- 
ficient). The term — 'jA 2 c in the right hand part of the Cahn-Hilliard equation regularizes solutions 
and the existence theorem was proved for the initial-boundary value problem given smooth initial 
data Il26ll . The proof relies essentially on the existence of a Lyapunov functional F (11.251) . 

The time derivative of F in a domain V with a smooth boundary and without external fluxes 
(( J, n) = on the boundary, where n is the vector of outer normal to the boundary) is 

dF 



\id t cdx = / /idivJdx 
v Jv 



, (1-27) 

= - / (V//, J) dx = -D / (V/i) 2 dx < . 
Jv Jv 

If we correct the Cahn-Hilliard equation by the "Teorell" factor c then the dissipation inequality 
F < (11.271) persists: Due to the Teorell formula 

f df c (c] 

j = -mcV/i = -mcV v ' - 7 Ac 



dc 

dc ( * f d n c ) a N x 

— = — div J = m div I c grad I — ^Ac 

Here, m is the Einstein mobility. This equation should be called "The Cahn-Hilliard-Teorell" 
equation. In accordance with (11.281) . 



dF 

— — = / fid t cdx — I /idivJdx 
>v Jv 



dt 



, (1-29) 

/ (V/i, J) dx = -m / c(V/i) 2 dx < . 
Jv Jv 



This dissipation inequality allows us to transfer all the results about solutions of the Cahn-Hilliard 
equation to the Cahn-Hilliard-Teorell equation. 



1.2.4. Teorell Formula for Non-perfect Systems 

It seems very natural that the flux is proportional to the concentration of particles: the average 
velocity is proportional to the force and the total flux is the product of the average velocity and the 
amount of moving particles. This could be proved in the framework of non-equilibrium thermo- 
dynamics and the theory of absolute reaction rates when the concentration of moving particles is 
small. In perfect gases or in dilute solutions the chemical potential is 

/i = RT In c + /io , 

where /i does not depend on c (it is a function of T and the state of the environment). In this 
case, we neglect the interaction between moving particles and use the Teorell formula (exactly as 
Einstein did in his theory of Brownian motion 30 years before Teorell Il24l ). 
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When the concentration of moving particles c is not small enough then the formula for perfect 
chemical potential is no longer valid and in front of c in the flux a special activity coefficient a 
appears. 

Such coefficients were introduced for diffusion by Eyring at al in 1941 ll66l and were used 
systematically for the theory of nonlinear diffusion by Gorban at al in the 1980-1986 [|46ll . Roughly 
speaking, the activity 



should substitute the concentration c in the Teorell formula with the proper renormalization of the 
mobility coefficients: 



The renormalized coefficient m' is defined by the condition: (m'a) / (mc) — > 1 for c — > 0, Vc — >• 0. 
This means that 



or we can write the Teorell formula for non-perfect systems using the usual Einstein mobility m 
defined for small concentrations and this standard value of chemical potential, /i : 



This formula is the main analogue of Fick's law for monomolecular diffusion in non-perfect media. 

For the Cahn-Hilliard-Teorell equation (11.281) . the Teorell formula for non-perfect systems 
(11.301) significantly changes the diffusion coefficient: the regularizing gradient term should appear 
in the activity coefficient. 

More details about activity coefficients in thermodynamics can be found, for example, in Chap- 
ter 9 of the classical book EE]. 

In the phase separation problem, the components are definitely non-perfect and the further 
correction of the Cahn-Hilliard-Teorell equation by the activity coefficients is necessary. 

The problem of the extension of the Cahn-Hilliard approach to multicomponent diffusion was 
discussed by various authors [1761 0. Elastic forces and plasticity are also taken into account 
[|4l|35l. Nevertheless, the problem of the proper equations of multicomponent nonlinear diffusion 
in highly non-homogeneous condensed phases is still open. From our point of view, there is no 
single "proper model" and the variety of possible models is very rich. In our work, we attempt 
to formulate the proper language for the description of the universe of these models similarly to 
chemical kinetics models. 

1.3. Main Ideas 

1.3.1. Mechanisms as Collections of Elementary Processes 

A complex process can be disassembled into several elementary processes. The dependence of 
the process rate (the flux) on the state (concentrations, chemical potentials and their gradients) 




J = m'a(— V/i + (external force per gram particle)) . 




J = m exp 




(—V/i + (external force per gram particle)) . (1.30) 
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is simple for elementary processes. The model of the whole process is assembled from these 
elementary "details". 

This idea was developed in chemical kinetics. In 1862-1867, Guldberg and Waage proposed 
the mass action law for equilibrium. In 1879 they developed the mass action law for dynam- 
ics. This idea was developed further by many researchers and after several dozen years it was 
transformed into a technology for the representation of complex processes: A complex reaction 
is represented as an ensemble of elementary reactions. The reaction rate has a simple monomial 
dependence on concentration. 

Van't Hoff HI 071 called the reactions that satisfy the mass action law "normal transformations" 
and found that "normal transformations take place very rarely". Now we can say that most reac- 
tions are complex and the interaction of several elementary reactions causes non-trivial complex 
("abnormal") behavior. 

Van't Hoff did not study complex reactions by disassembling them into several elementary 
reactions. Therefore, he was disappointed with the mass action law and finally wrote: "As a 
theoretical foundation I did not accept the concept of mass action, I had to abandon this concept in 
the course of my experiments...". ("J'ai adopte pour la theorie, non la notion des masses actives, 
notion que j'ai du abandonner dans le cours de mes experiences, ..." H107H , P- 7.) 

The set of elementary reactions which constitute a complex reactions is called the reaction 
mechanism. A mechanic analogue is obvious: the elementary reactions are the details of the mech- 
anism that represents the complex reaction. This notion was, finally, introduced into chemistry in 
the 20th century due to the efforts of M. Bodenstein. 

M. Boudart described the "century of Bodenstein" in his paper (5J: "First came the data, then 
the rate equation, and finally the fitting of the data into the rate equation by means of a hypothesized 
mechanism with rate constants chosen for the best fit." 

The great success of this approach was the theory of chain reactions. N.N. Semenov was 
awarded the Nobel prize for this theory [|94| . The modern theory of complex chemical reactions 
is based on the idea of the detailed reaction mechanism and the simple kinetic law of elementary 
reactions H1171I . 

Many authors proposed various particular mechanisms of nonlinear diffusion. One of our goals 
in this work is to repeat the way of chemical kinetics in application to multicomponent diffusion 
and to create a comprehensive theory of the mechanisms of diffusion. 

1.3.2. Discrete Kinetic Models and Lattice Automata 

In the 1940s, S. Ulam and J. von Neumann proposed networks of interconnected finite-state au- 
tomata for the modeling of complex systems. In the first period of study, research was focused on 
the abilities of these networks and, in particular, on the ability of self-reproduction HI 111 . 

The behavior of these cellular automata is so variable and surprising and its complexity is so 
high that Ulam proposed the idea of the computational experiment: we should regard our invention 
as a new sort of reality and study it by the experimental approach, as physics or chemistry does. 

Cellular automata were invented as an intellectual journey but soon were recognized as efficient 
tools for modeling [110311 . 
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Feynman's attention to automata with local interactions as a tool for simulating physics, at- 
tracted much attention to this area: "Therefore my question is, can physics be simulated by a 
universal computer? I would like to have the elements of this computer locally interconnected, and 
therefore sort of think about cellular automata as an example (but I don't want to force it). But 
I do want something involved with the locality of interaction. I would not like to think of a very 
enormous computer with arbitrary interconnections throughout the entire thing" OTTl . 

The idea of modeling the natural world in terms of the behavior of sets of rules that can be em- 
bodied in simple automata with local interactions is now an important part of science. Sometimes 
it is called "the new science" to distinguish this approach from classical modeling by equations 

ma. 

For the modeling of transport processes the lattice gas automata H114[ [T6l were invented and 
the lattice Boltzmann methods ||98l became very popular: they are flexible and efficient. At the 
same time, the lattice Boltzmann methods are very simple for programming and parallelization. 

The essence of the lattice Boltzmann methods was formulated by S. Succi in the following 
maxim: "Nonlinearity is local, non-locality is linear" H100II . We should even strengthen this state- 
ment. Non-locality (a) is linear; (b) is exactly and explicitly solvable for all time steps; (c) space 
discretization is an exact operation. 

The lattice Boltzmann method is a discrete velocity method. The finite set of velocity vectors 
{vi} (i = 1, ...m) is selected, and a fluid is described by associating, with each velocity Vi, a 
single-particle distribution function fa = fi(x,t) which is evolved by advection and interaction 
(collision) on a fixed computational lattice. The values ft are named populations. If we look at 
all lattice Boltzmann models, one finds that there are two steps: free flight for time 5t and a local 
collision operation. 

The free flight transformation for continuous space is 

fi(x, t + St) = fi(x - ViSt, t). 
After the free flight step the collision step follows: 

fiW^Fiiiftix)}), (1.31) 

or in the vector form 

f(x) F(f(x)). 

Here, the collision operator F is the set of functions Fj ({/_/}) (i = 1, ...m). Each function Fj 
depends on all fj (j = 1, ...m): new values of the populations fa at a point x are known functions 
of all previous population values at the same point. 

The lattice Boltzmann chain "free flight — > collision — > free flight — > collision • • • " can be 
exactly restricted onto any space lattice which is invariant with respect to space shifts of the vectors 
ViSt (i = 1, . . . , m). Indeed, free flight transforms the population values at sites of the lattice into 
the population values at sites of the same lattice. The collision operator (11.311) acts pointwise at 
each lattice site separately. Much effort has been applied to answer the questions: "how does the 
lattice Boltzmann chain approximate the transport equation for the moments Ml", and "how does 
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one construct the lattice Boltzmann model for a given macroscopic transport phenomenon?" (a 
review is presented in the book (981). 

The lattice Boltzmann models should describe the macroscopic dynamic, i.e., the dynamic of 
macroscopic variables. The macroscopic variables Mt(x) are some linear functions of the popu- 
lation values at the same point: Mi(x) = ^ maf^x), or in the vector form, M(x) = m(f(x)). 
The macroscopic variables are invariants of collisions: 

X>*/< = Y^rnaFidfj}) (orm(/) = m(F(f))). 

i i 

The standard example of the macroscopic variables are hydrodynamic fields (density-velocity- 
energy density): {n,u, E}(x) := ^{1, vf/2}fi(x). But this is not an obligatory choice. On 
the other hand, the athermal lattice Boltzmann models with a shortened list of macroscopic vari- 
ables {n, nu} are very popular. 

The quasiequilibrium is the positive fixed point of the collision operator for given macroscopic 
variables M. We assume that this point exists, is unique and depends smoothly on M. For the 
quasiequilibrium population vector for given M we use the notation f^, or simply /*, if the 
correspondent value of M is obvious. We use II* to denote the equilibration projection operation 
of a distribution / into the corresponding quasiequilibrium state: 

n*(/) = /; (/) . 

Usually, collision operators are taken in the form: 

F(f):=U*(f)+A(U*(f)-f), (1.32) 

where A is a linear operator, whose spectrum belongs to the interior of the unit circle. A special 
case of (11.321) is very popular, the lattice Bhatnagar-Gross-Krook (LBGK) model: 

F(f):=f + u(U*(f)-f). (1.33) 

In this brief introduction of LBM we follow the paper §8§. 

The simplest LBGK realization of Fick's law (in ID) gives the following system. The discrete 
velocity set includes two elements only, v and —v. The time step is r the correspondent grid step 
is h — vr. The microscopic variables, the populations, are: /~ for velocity — v and / + for v. The 
macroscopic variable, the density, is p = f~ + / + . The corresponding equilibrium is U*(f) = /*: 

f~ + f + 



For the non-negative populations, the equilibrium distribution is the maximizer of the entropy 
S = —f~ In /~ + / + In / + under a given value of the macroscopic variable p. 
Let us take the LBGK collisions (11.331) with to = 1, i.e. 

F(/) ± = n*(/) ± = ; (1-34) 
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This particular case of the LBGK collision integral is an example of the so-called Ehrenfests' 
coarse-graining. The idea of artificial partial equilibration steps was proposed by T and P. Ehren- 
fest for the foundation of statistical physics [23] and further developed to a general formalism of 
nonequlibrium thermodynamics ll48l[T8ll4~9ll . 

A review and comparative analysis of different approaches to coarse-graining was published in 

ma. 

The LBGK chain for the collision integral (11.341) has a very simple form: 

f + (nh, (m + l)r) = /"(n/i, (m + l)r) 
}' + {{n — l)h, mr) + f~((n + l)h, mr) 

(1.35) 

~ 2 ' 

Therefore, for the density p we get 

P (nh, (m + l)r) = P(( " " 1)h > mr) ± P(( " + 1)h > mr) . (1.36) 

This appears to be the one of the most common explicit finite difference methods for Fick's diffu- 
sion equation. The diffusion coefficient is D = h 2 /(2r) = v 2 t/2 and depends explicitly on the 
lattice parameters. We can decouple D and the lattice parameters if we use u E [1, 2] in the LBGK 
collision integral (11.331) . This lattice-gas scheme does not coincide with any of the finite difference 
schemes. Nevertheless, it also models diffusion and, to the first order in r, D = v 2 r^f- ||98ll44l . 

Now, the area of applications of the cellular and lattice Boltzmann automata is very wide and, 
in addition to classical fluid dynamics, includes many areas of chemistry [|65l . models of phase 
separation ll92l . dynamics of macromolecules and many other topics. 

We use cellular automata and lattice models of nonlinear multicomponent diffusion for two 
purposes: 

• As a tool for model creation (after that, this model could be translated into other languages, 
such as partial differential equations (PDE); 

• As a tool for numerical simulation without the stage of PDE model. 

Elliott and Stuart lETll used the cell model of diffusion to study semilinear parabolic equations. 
They proved the existence of absorbing sets, bounded independently of the mesh size for discrete 
models. Discrete Lyapunov functions were constructed. We use the special quasichemical ap- 
proach for the generation of the the cell models 11461 that allowed us to construct the Lyapunov 
functions for semi-discrete systems and to prove stabilization of the solution in space and time 
under proper conditions. 

Let us consider our space divided into cells, a system represented as a chain of cells of ho- 
mogeneous composition and elementary transfer processes between them. It is sufficient for our 
purposes to discuss two sells (Fig [2]). Let us numerate these cells by the Roman numbers I and 
II and mark all the components and quantities related to them by the upper index I or II, corre- 
spondingly. The lists of the components for cells are different just by the upper index: A\, . . . A T n , 
A 11 A 11 
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Figure 2: Cell Jump Model 



The mechanism of diffusion is defined as a list of elementary transitions between cells de- 
scribed by their stoichiometric equation. Since diffusion is a sort of jumping reaction on the border, 
for these jumps the stoichiometric equation is written as, 



E ««4 + E E + E ^ 



(1.37) 



where r is the number of processes, a*'/ 1 , and /3^ n , are the stoichiometric coefficients which in- 
dicate the number of particles in cells involved in the process. The direction of changes in the 
elementary event (11.371) is defined by two stoichiometric vectors 



7r 



. ' r% 5 in 



8 n - a 11 



Examples of elementary acts are presented in Fig.[3l 

Elementary events (11.371) should not include reactions. Therefore, for each i, the amount of A^ 
in the system (A\ + A] 1 ) should not change. This means exactly that for all z, r 



-7° 

in 



Let us use the notation 



It 



-7 11 

Iri 



The composition of each cell is a vector N 1 ' 11 . The components of this vector, Nj' u are the 
amounts of Ai in the correspondent cell. We describe the dynamics of the compositions of two 
cells by the equations: 

dN l dN* ^ 

b lr w ri (1.38) 



dt 



dt 



where S is the area of the boundary between two cells and w r is the rate of the process. For many 
cells the equations are the same, but with more pairs of cells interacting, and therefore there are 
more terms. 

The rates are intensive variables and should be defined as functions of concentrations or chem- 
ical potentials. The crucial question is: how to describe function w r (c l , c 11 ), where c I,n are con- 
centrations components in cells. 

The real physics of diffusion may be more complicated. For example, the intensity of jumps 
and the reaction rate w r (c I , c 11 ) may depend not only on (c 1 , c 11 ) but on the surrounding. For 
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(a) Simple diffusion: a particle from the cell I jumps into the cell I] 
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(b) Jumps to free places: a 
place in cell II and inverse 
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(c) Jumps with clustering: two particle attract the third one 



Figure 3: Elementary acts of diffusion, examples. 
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Figure 4: Cell Jump Model with first surroundings. 

example, direct simulation of the jumps on the surface flU demonstrates that the influence of the 
surrounding is crucial for structures and critical effects on the surface. 

For each process (11.371) there is the space-inverted process that is defined just by changing I 
to II and vice versa. We mark the quantities for the space-inverted processes by '. For example, 
7' = —7. The detailed space-inversion symmetry requires that the rate functions for them should 
differ just by the transposition of the vectors of variables, c 1 , c 11 : 

w' r (c 1 ,c 11 ) = w r (c I1 ,c 1 ). (1.39) 

This requirement of detailed space symmetry allows us, in particular, to exclude various types of 
advection and transport driven by external force. Diffusion, by its definition, is driven by the 
gradients of the concentrations (or, in the thermodynamical approach, by the gradients of the 
chemical potentials). This is not the only way to formulate of pure diffusion equations without 
advection. Another possibility gives us, for example, the diffusion systems with complex balance 
(Section [23]). 

There are three ways to define the rate functions: from a phenomenological law (like the mass 
action law), from thermodynamics (like the generalized mass action law) or by direct stochastic 
simulation of particles jumps in cells (like in the Gillespie approach [|40ll4T| ). 

In our research, we focus on the first two approaches. Therefore, we consider our lattice model 
as a semi-discrete model (discrete in space and continuous in time). For this semi-discrete model, 
the system of kinetic equations (11.381) describes diffusion. The continuous limit of these equations 
gives us the diffusion PDE. The discrete scheme by itself can serve as a computational model. 

A couple of simple examples can clarify our approach: 



Simple diffusion (Fig.|3(a)j), A\ — > Af and Af — > A\ with the same rate constants. Particles 



jump into the neighbor cells. For perfect mixtures, w r = kc\, w' r = kc\ l and in the contin- 
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uous limit we get Fick's law (11.2b as the first Taylor approximation. In this approximation, 
D = kl where I is the cell size. 

• Jumps to free places (Fig. [3(b)]), A\ + Z u -> Af + Z 1 and Af + Z 1 ->■ A\ + Z u . According 
to the mass action law, w r ((?, c 11 ) = kc\z u , w'^c 1 , c 11 ) = kcfz 1 , where z is the concentration 
of free places. In the first Taylor approximation, J = —kl(zVci — CjVz) and we get the 
model (071) . (fU9l) . 

To get the continuous limit, we take c 1 = c(x), c 11 = c(x + I) and use the Taylor expansion: 
c(x + I) = c(x) + ld x c + o(l). If we the consider a sequence of cell representations of diffusion 
with various I then, for the invariance of the first order, the scaling rule should be implemented: 
D = kl does not change with the change of size, therefore, the rate constant k depends on I: 
k = D/l. 

It is not always possible to keep to first order only. If this approach gives a negative diffusion 
coefficient then for regularity we have to keep the higher derivatives. For example, let us take the 
diffusion mechanism with attraction: 

A] + 2Af ^3Af. (1.40) 

The space-inverted process in this case does not coincide with the inverse one. If we change the 
upper indexes (I to II and II to I) then we obtain 

2A\ + A? ^ 3A\ . (1.41) 

This mechanism means that 2 particles attract the third one. This mechanism is represented in 



Fig. 3(c) 



The reaction rates are: 



w r = k r c\{cY)\ w' r = k r {c l ^c l l 



The flux of Ai from the first cell to the second one is 



J = W r — w' r = krclcf^f — C-) 



Therefore, to first order we have 



J = klc 2 Vc = -klVc 3 ; 
3 



the sign is opposite to standard diffusion. This flux goes in the direction of gradients. The diffusion 
equation is 

8c 1 

— = -fc/div(c 2 Vc) = -kl-Ac 3 . (1.42) 
Of course, if we take the mechanism (n > 1) 

A] + nAf ^(n + l)A\ , A\ + Af -> (n + 1)^ , 
then we get the equation 
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CfC T) — 1 

^ = -kl(n - l)div(c n Vc) = -kl Ac n+1 . 

at n + 1 

This diffusion process has two properties: first, it goes along gradients and all deviations from 
the uniform state will increase. Second, this diffusion is slow for small concentrations (the diffu- 
sion coefficient goes to when c approaches 0) and accelerates with the concentration growth. 

The equation d t c = —DAc n (n > 1) admits a family of self-similar solutions with bounded 
support, which collapse in finite time. These solutions have the form 

A , (r 



C(T) = —0 1- 



where 

• r is the time till collapse; 

• q is the dimension of space (usually, q =1, 2 or 3); 

• p is the radius of the sphere, outside of which the solution is zero 

p = B(Dt) 8(«-!)+2 ; 



• = (1 -§ 2 )^ for?? < 1 and0($) = Oiftf > 1; 

• The constants A, B depend on q, n and the total amount N = J c(x) dx. 

This is the so-called Barenblatt solution Q]| for the equation of porous media d T c = +DAc n . Such 
solutions were used in the analysis of an explosion which starts from a singularity for equations 
d t c = +DAc n (the classical review of self-similar solutions was published by Barenblatt and 
Zeldovich [H). 

The cell model of diffusion with attraction (11.401) for a finite number of cells of a given size / 
is a rather regular system of nonlinear ODE, but to first order of the Taylor expansion in / the PDE 
(11.421) produces a singularity in an arbitrarily short time from smooth initial data. The second order 
Taylor approximation adds nothing because the even terms in / cancel out if we take into account 
both the left and right neighbors of the cell. The third order Taylor expansion gives a regularized 
equation: 

d ( I 2 d 2 c\ 
J = J = Wr - w i=klc*-{c + -—yo{?)- 

d ± -kl—c 2 — (c 

dt dx dx \ 3 dx 2 J 

This is an example of the Cahn-Hilliard type equation for spinodal decomposition with the regu- 
larizing term — div(c 2 gradAc). In this equation, the cell size cannot be eliminated by scaling. The 
length / is the "regularization length". All inhomogeneities of size smaller than / are smoothed by 
the biharmonic term. 
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As we can see, the mass action law and the cell representation of the elementary acts of dif- 
fusion give the opportunity to model the Cahn-Hilliard type phase separation. Nevertheless, the 
approach based on the non-perfect thermodynamic potential (11.251 ) gives a better representation of 
the basic physics and does not require complicated elementary processes. Just the simplest Fick 
scheme, 

A-i — >• A i , A i — > A i 

with the non-perfect Ginzburg-Landau free energy gives the Cahn-Hilliard equation (Sec. 13.2.1) . 
The diffusion mechanism with attraction (11.401 ) (Fig. |3(c)[ ) differs from the elementary Fick 



mechanism (Fig. |3(c)| ) and from the mechanism of jumps to free places (Fig. |3(b)] ). The dynamical 



difference is obvious, the diffusion mechanism with attraction generates instabilities of the homo- 
geneous state, clustering and singularities. On the other hand, Fick's law and the mechanism of 
jumps to free places (Fig. |3(b)| ) allow a global Lyapunov functional and, in the systems without 
external fluxes, lead to homogeneous equilibrium. 

These mechanisms have also a very important structural difference. If we look at the direct 
and the space-inverted processes (Fig. [3]) then we find that for the first two mechanisms, the space- 
inverted processes coincide with the inverse processes, which we get just by inversion of the arrow 
(or by the exchange a and (3 coefficients in the stoichiometric equations (11.37I )). For the elementary 
processes with attractions (Fig. |3(c)[ ) the inverse processes are processes with repulsion: 



3 A? ^A\ + 2Af , 3 A] -> 2 A] + Af . (1 .43) 

The diffusion processes for which space-inverted elementary processes coincide with the in- 
verse processes, have a fundamental property: The entropy production is positive for the corre- 
sponding mass action law diffusion equations. 

Let us consider a complex diffusion process in a bounded domain with smooth boundary and 
without external fluxes. We prove the following theorem in Section |2.2.3.| 

Theorem 2. Let a complex diffusion process consist of elementary processes, which satisfy 
the following property: the space-inverted elementary process coincides with the inverse process. 
Then, for the mass action law equation of diffusion ( 12.251 ), the principle of detailed balance is 
valid, the global convex Lyapunov functional exists and the uniform distribution is asymptotically 
stable. 

This global Lyapunov functional may be selected in the form of the (minus) classical entropy, 
the sum of terms c In c for all cells and components, or, for the continuous limit, 



d In Cj dx 



A particular case of the dissipation inequality for such processes is inequality (11.221) for the 
diffusion equations (11.201) that describe diffusion by exchange of positions. It is valid because the 
exchange mechanism satisfies this fundamental property: the space-inverted elementary processes 
coincide with the inverse processes. 
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1.3.3. Thermodynamics and Intermediate Complexes 

Thermodynamics is not always a good leader, but it is always a good judge. We cannot create 
nonlinear equations directly from thermodynamic principles, but we must always check whether 
our equations satisfy thermodynamics. They should satisfy the thermodynamic restrictions if we 
do not want to produce a perpetuum mobile in our theory. 

We also include some other fundamental restrictions like micro-reversibility in the thermody- 
namic requirements. 

It is not always simple to coordinate the lattice models with thermodynamics, nevertheless it is 
possible 116411991 . 

There are two main approaches for the introduction of thermodynamics into kinetic models. 
First, we can start from general kinetic equations based on the representation of a complex process 
as an ensemble of elementary processes with a given simple kinetic law of elementary processes 
(for example, the mass action law). After that, we will find that the rate constants of the ele- 
mentary processes are not independent. They must be coordinated to meet the thermodynamic 
requirements. Therefore, not all the possible kinetic systems are allowed thermodynamically. 

Another approach starts from the thermodynamic description of the system. We should find 
thermodynamic potentials which describe the system under given conditions. We have to know 
entropy, free energy (Helmholtz energy), or free enthalpy (Gibbs energy) for the proper set of 
independent variables lfT4l . After that, we define the rate of elementary process through the ther- 
modynamic functions but with some arbitrariness: some constants remain free of thermodynamic 
restrictions. These constants are independent for different elementary processes. 

Which way is better? It is not a proper question: both are good for their purposes. The first 
approach (we start from kinetics and then add thermodynamics) is very flexible. In particular, it 
can be used when thermodynamic restrictions are not needed. For example, when we consider sub- 
systems of open systems like the system of surface components in heterogeneous catalysis, then 
the constants of elementary processes include additional dependencies on some additional concen- 
trations and are not the "proper" rate constants. Therefore, they do not satisfy the thermodynamic 
restrictions, and a subsystem may demonstrate non-thermodynamic behavior like non-decaying 
oscillations or bifurcations. 

The second approach is unavoidable for non-perfect systems. The kinetic law of elementary 
processes depends on the thermodynamic potential. For all perfect systems it is the same mass 
action law, but any deviation from the perfect thermodynamic function requires its own deviation 
of the kinetic law from the mass action law ll68l . This deviation may be reformulated as the 
generalized mass action law with activities instead of concentrations but the activities are defined 
through thermodynamic potentials. 

In our work we follow both approaches: First, we formulate the mass action law for diffu- 
sion and study this with and without the thermodynamic restrictions. Secondly, we introduce the 
thermodynamic formalism for diffusion in non-perfect systems. The ideas for both approaches for 
diffusion were formulated in the early 1980s Il53ll46l . The detailed analysis of the thermodynamic 
restrictions on chemical kinetics was performed by Gorban in 1982 Il42l . 

Our work was influenced by the works of N.G Van Kampen H105H , M. Feinberg [|29l[30l and 



24 



A.N. Gorban et al. 



Quasichemical Models of Multicomponent Diffusion 



Horn and Jackson ||62l . In 1973, N.G. Van Kampen proposed a general formulation for the rates of 
irreversible processes as a combination of "unilateral transfer flows". Each unilateral flow transfers 
energy and particles in one direction. Van Kampen decomposed the total in partial systems, each 
of which is in equilibrium and therefore possesses a well-defined temperature, entropy, and other 
thermodynamic quantities. Although the total system Y is not in equilibrium, it is still possible to 
attribute an entropy to it. Then Van Kampen studied the unilateral fluxes between subsystems. 

We decompose the Van Kampen unilateral processes further and represent them as a collection 
of essentially one-dimensional elementary processes with the simple kinetic mechanism, the mass 
action law or the generalized mass action law. 

We start from a similar representation of the total system and supplement it with the system 
of stoichiometric equations of elementary unilateral processes. To find the rate of the elementary 
processes we use an idea of intermediate complex (compound). This approach is borrowed from 
the theory of absolute reaction rates but we do not use the special idealization of the reaction pass 
and postulate the more general microscopic Markov kinetics instead. 

If the concentrations of compounds are small and the equilibrium between intermediates and 
other components is fast (both assumptions are important) then we approach the generalized mass 
action law, which is very similar to the Marselin-de Donder kinetics and the generalized mass 
action law studied by Feinberg, Horm and Jackson and other authors ||29l |3~0[|62[ [TOl . 

This formalism is very convenient for implementation of the microreversibility consequences 
in the form of detailed balance conditions 11731 . In addition, if there is no microreversibility then 
the thermodynamic behavior is also guaranteed by the special more general relations between ki- 
netic constants, which follow from the Markov kinetics of intermediate complexes. First, the idea 
of such relations was proposed by Boltzmann as an answer to the Lorentz objections against Boltz- 
mann's proof of the //-theorem. Lorentz stated nonexistence of inverse collisions for polyatomic 
molecules. Boltzmann did not object to this argument but proposed the "cyclic balance" condi- 
tion, that means balancing in cycles of transitions between states Si — > S2 — > ■ ■ ■ — > S n — > Si. 
Almost 100 years later, Cercignani and Lampis [15] demonstrated that the Lorenz arguments are 
wrong and this Boltzmann new relations are not needed for the polyatomic molecules under the 
microreversibility conditions. The detailed balance conditions should hold. 

Nevertheless, this Boltzmann's idea is very seminal. It was studied further by Heitler ll60l 
and Coester |[T9l and the results are sometimes cited as the "Heitler-Coestler theorem of semi- 
detailed balance". In 1952 [97] proved these conditions for the Boltzmann equation. For the 
micro-description he used the S'-matrix representation, which is in this case equivalent for the 
Markov microkinetics (see also H112II ). Later, this sort of relation was rediscovered for chemical 
kinetics Il3~0~ll62~l . The general proof for nonlinear nonequlibrium processes was presented recently 
ll54l . In our analysis of these Boltzmann-... -Stueckelberg relations we follow the later. 

We extend the usual stoichiometric equations by additional reactions: an input linear com- 
bination of reagents forms a corresponding compound; this compound transforms into another 
compound that disintegrates into the corresponding output linear combination of reagents: 

i i 

Here p is the elementary reaction number. 
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Figure 5: A 2n-tail scheme of an extended elementary process (11.441) . 

It is useful to visualize the reaction scheme. In Fig. [5] we represent the 2n-tail scheme of an 
elementary reaction sequence (11.441) . This scheme was proposed in [|54l . 

We assume that the amount of each compound B p is small enough to apply the perfect entropy 
formula, and that the equilibrium between each compound and the corresponding linear combi- 
nations of reagents is fast enough to apply the quasiequilibrium approximation [|54l (the detailed 
analysis of this approximation was given in |j52l[50l[5T1l ). 

The main difference from the Eyring approach ll66l is in the model of the hidden reaction of 
the "activated complexes": 

• Eyring used for each reaction one complex with a continuum of energetic states along the 
"reaction path", whereas we use two compounds (or two states); 

• Eyring modeled the reaction of the intermediate complex as a classical motion along the 
additional coordinate and even added this degree of freedom with classical kinetic energy of 
this motion to the free energy calculation, whereas we follow the Stueckelberg approach and 
model this reaction as a first order Markov kinetics, a Markov process with two states. 

The difference in the macroscopic consequences of these approaches seems to be not very large 
because the result of the Eyring approach is one relaxation time approximation for each reaction. 
From the dynamical point of view, this result coincides with the two-state Markov model. 

The main differences may arise in the hints which these approaches give to the microscopic 
calculation of the macroscopic quantities. In our work, we concentrate on the macroscopic dy- 
namics. 

2. Mass Action Law for Diffusion 

2.1. Mass Action Law 

2.1.1. Mass Action Law Kinetic Equations 

This is an auxiliary subsection where we collect main definitions and results about the Mass Action 
Law (MAL). 
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To construct a system of kinetic equations by MAL, one needs the following inputs: 

1 . A list of components; 

2. A list of elementary reactions represented by their stoichiometric equations; 

3. A set of reaction rate constants. 

The list of components is just a set of symbols (the component names). We usually assume that 
this set is finite, Ai, A 2 , . . . , A n . 

Elementary reactions are given by their stoichiometric equations, 

J2 a nAi^J2PriAi, (2.1) 

i i 

where r is a reaction number, a ri and ri are nonnegative numbers, the stoichiometric coefficients. 
By default, they are assumed to be integer but, sometimes, there occurs a need in nonnegative real 
coefficients. 

For each elementary reaction (12.11) . a stoichiometric vector is defined, 

TV • Tri 0ri O^ri • 

This is a "bookkeeping" vector, whose components are "gain minus loss" (or "income minus out- 
come"). 

We will also use the loss and gain vectors of elementary reactions: a r (loss) with coordinates 
a ri and /3 r (gain) with coordinates /3 ri . Of course, 

7r = Pr - 0> r . 

The stoichiometric matrix T is the matrix with columns j^. = 7^, the first index in 
corresponds to component and the second index corresponds to reaction. Reaction rate constants 
k r are non-negative numbers. They should be defined for all elementary reactions. For each com- 
ponent Ai, a real variable, concentration q is defined. Vector of concentrations c has coordinates 

c { . 

The reaction rate for the elementary reaction (12.11) is the following function of c 

n 

W r = k r Y[ . (2.2) 

i=l 

The MAL kinetic equations are 

— = ^2 ^rWr ■ (2.3) 

r 

From the physical point of view, these equations describe isochoric isothermal processes for 
perfect systems. For non-isochoric or non-isothermal processes it is necessary to introduce the vol- 
ume (together with pressure) and the enthalpy (together with temperature) explicitly and describe 
their dynamics. 
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For a given reaction mechanism, a linear stoichiometric conservation law is a linear functional 
b(c) = biCi that annihilates all stoichiometric vectors: 

Kir) = for all r . 

The stoichiometric conservation law is strictly positive, if all bi > 0. The assumption about exis- 
tence of a positive stoichiometric conservation law plays an important role in the MAL kinetics. 

2.1.2. Existence and Uniqueness of Solutions 

Let in the reaction mechanism all nonzero coordinates of the loss vectors be not less than 1: 

a ri > 1 or a ri = . 

This assumption is valid, for example, if all the stoichiometric coefficients are nonnegative integers. 
Let us assume also that there exists a strictly positive stoichiometric conservation law b. Then the 
following existence and uniqueness theorem for the MAL equation holds. 

Theorem 1. For any nonnegative initial data c(0) fcj(0) > 0) there exists a unique solution of 
A2.3\) c(t)for all t > 0. This solution is nonnegative (Ci(t) > 0) and satisfies the conservation law: 
b(c(t)) = 6(c(0)). 

This is a well known result (see, for example, H109II ). The proof is quite simple. First, of all, let 
us consider a bounded neighborhood U of the simplex S : c% > 0, b(c) = 6(c(0)). The right hand 
site of the MAL kinetic equations (12.31) has continuous first derivatives and these derivatives are 
bounded in U . Therefore, according to a standard existence and uniqueness theorem its solution 
exists on some time interval t E [0, T] and this T is the same for a compact set of initial data 
c(0) 6 S d [/. Secondly, let us mention that if j ri < then a ri > 1 and the reaction rate w r (12.21 ) 
includes the factor c"" . Therefore, 

^2 IriWr = Cig(c) , 

7ri<0 

where g(c) is continuous function. 

If Cj — > then X] 7ri <o 7« w r — > 0. If q = then q = Y^rlriW,- > 0. Therefore, the simplex 
S is positively invariant with respect to equations (12.31 ): the existent solutions c(t) do not leave 
this simplex for t > if c(0) e E . Finally, this implies global existence of solutions in E . 

Both conditions of existence of a strictly positive stoichiometric conservation law b and of 
coordinates of the loss vectors, «j = or oti > 1 are significant. Just for example, we can consider 
mechanisms that violate these conditions: 2A — > 3A and \A — > A. For the first mechanism, 
c = kc 2 and there is no global existence, for the second system, c = ky/c and there is no uniqueness 
of solution. 
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2.1.3. Detailed Balance 

The expected behavior of a system of physical or chemical kinetics is simple in the absence of 
external fluxes: everything goes to equilibrium. A Lyapunov function for this relaxation is the 
corresponding thermodynamic potential. 

The MAL kinetics (12.31 ) do not assume any thermodynamic properties "from scratch". More- 
over, this class of kinetic equations is so rich that it is dense in the class of all smooth semi- 
dynamical systems in S for a given conservation law b [46]. Additional assumptions are needed 
to guarantee the thermodynamic behavior. 

The most celebrated sufficient condition for the thermodynamic behavior of the MAL kinetics 
is the principle of detailed balance. This principle, as a realization of micro reversibility was known 
for the Boltzmann equation |[6l (since his proof of the //-theorem in 1872) long before Onsager's 
reciprocal relations If83ll84l . A. Einstein used this principle for the linear kinetics of emission and 
absorption of radiation ll2~5i In 1901, R. Wegscheider published analysis of detailed balance for 
chemical kinetics [11 1311 . 

To formulate the principle of detailed balance, it is convenient to join pairs of direct and inverse 
elementary reactions in (12.11) and write 

ariAi ^ Yl P"Ai . (2.4) 

i i 

If the inverse reaction does not exist in the original mechanism, we formally add it but assume 
that its rate constant is zero. We mark quantities for the direct and inverse reactions by the upper 
indexes + and ~ and write the MAL reaction rate: 

W r = — W~ , 

w+ = k+f[dT, w; = k;f[ct. (2 ' 5) 

i=l i=l 

For the MAL kinetics the principle of detailed balance is: there exists a strictly positive point 
of detailed balance that is such vector of concentration c* that c* > and 

w+(c*) = w;(c*) (= w* r > 0) . (2.6) 

This means that at least at one positive point the direct elementary processes are equilibrated by 
the inverse elementary processes. 

Existence of one such point implies that all equilibria are also points of detailed balance (12.61) 
and, moreover, there exists a global Lyapunov function that has the form of relative entropy. This 
is, precisely, the analogue of the Boltzmann //-theorem for the MAL kinetics. 

For the formulation, use and proof of this theorem, it is convenient to rewrite the formulas for 
the direct and inverse reaction rates (12.51 ) using c*, w* and the detailed balance relations (12.61) : 

W r = Wr — W~ , 
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The Lyapunov function is: 



G = E^( ln (|)- 1 )+E c *- ( 2 - 8 ) 



Here, the last constant term is added to satisfy G(c*) = 0. 

The partial derivatives of G (the analogs of chemical potentials) are 



dG , / a , 

i-=ln 4 ■ (2-9) 



dot \ c 

Therefore, we have one more form for the MAL kinetic law with detailed balance (12.71) : 



w*exp J ^°Vi^— ) = w* exp(ov, V C G) 



<9q / (2.10) 



5^ &r ) = ^ exp ^ r ' 



It is worth mentioning that 



to 



^ = exp[-( 7r ,V c G)]. 



For the time derivative of G due to the MAL kinetics (12.31) with the detailed balance, simple 
algebra gives the dissipation inequality: 



dG 



J2w r ( lr ,V c G)^-J2(^r-w;)hi(^j <0. (2.11) 



The last inequality holds because In x is a strictly monotone function and lnx — In y has the same 
sign as x—y has. Obviously, G\ c = if and only if c is a point of detailed balance. This equilibrium 
point may be different from the point c*, which was used for the definition. All the positive points 
of detailed balance for the MAL system (12.71 ) form a smooth manifold with dimension 

n-rank{7 1 ,7 2 ,...}, 

where n is the number of components and rank{7i, 72, . . .} is the rank of the system of the stoi- 
chiometric vectors for the given reaction mechanism. 

If we fix values of all stoichiometric linear conservation laws then the strictly positive point 
of detailed balance is unique for the given values. Indeed, the dissipation inequality is valid for 
every pair of mutually inverse reactions and all the terms uv( 7r , V C GQ in G (12.1 II) are non-positive. 
Therefore, at any strictly positive equilibrium point c cq , (7,., V C G) = for all r. This means that 
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c eq is a critical point of G in c eq + span{7i,7 2 , . . .}. The function G is strictly convex at any 
positive point: its Hessian is positive definite, 

d 2 G 1 

dcidcj q lJ ' 

where Sy is the Kronecker delta. Therefore, there may exist only one positive critical point of G 
on a linear manifold. 

Due to the logarithmic singularity of the gradients at the boundary of M + (where some of 
Ci = 0), G achieves its global minimum in 

(c eq + span{ 7 i,72,...})n M + 

at a positive point. This point is a positive point of detailed balance. 
For any positive vector, c°, the polyhedron 

V = (c + sp a n{ 7l ,7 2 ,...})f|R';, 

is positively invariant with respect to (12.3b . For systems with detailed balance it includes one and 
only one positive point of detailed balance. This was first demonstrated by Zeldovich in 1938 
(reprinted in 1996 H118I0 . In our analysis we mainly follow II109II . 



2.1.4. Complex Balance 

In this subsection we consider the direct and inverse reactions separately as we did it before in 
(12.11) . (12.21) . Detailed balance is s sufficient but not necessary condition of the thermodynamic 
behavior. The simple example of thermodynamic behavior gives any monomolecular (linear) re- 
action mechanism, which consists of reactions A, — > Aj. Let us use notation kji for this reaction 
rate constant. 

The MAL equations for a monomolecular reaction mechanism are 

—y- = (hjCj — kjiCi) . (2.12) 

Let c° be a strictly positive steady state for these equations (not necessarily a point of detailed 
balance): 

With this c° we can rewrite the second term in (12.121) : 

c° 

kjj = k*i^ ; 
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Therefore, the kinetic equations (12.121) have the equivalent form for given c°: 



dt 



Then we can define 



After simple transformations, we find that due to (12.131) 

dG 



iS(S)-0-|(H| 



< o. 



To prove this formula, it is worth mentioning that for any n numbers ft, 

E ^(/*-/i) = o- 

This gives us the first two terms in the square brackets with 



The last term, 



cyy \cj cy 



(2.15) 



appears in the straightforward computation of the time derivative of G due to kinetic equations 
(I213T) . 

The expressions in square brackets in (12.151) have the form 

f(a)-f(b) + f'(a)(b-a) 

for the convex function f(x) = x(\nx — 1). This expression is always non-positive because of 
Jensen's inequality. 

Linear MAL kinetics can obviously violate the principle of detailed balance. For example, an 
irreversible cycle 

Ai A 2 -> . . . ->■ A n -> Ai 
has always a positive steady state but never has a positive point of detailed balance if n > 2. 
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There is a nice generalization of the dissipation inequality ( 12.151) for the nonlinear MAL equa- 
tions under some algebraic conditions on the kinetic constants. These conditions are strictly 
weaker than the principle of detailed balance. They were discovered for the Boltzmann equation 
by Stueckelberg [97] in 1952 and called later the "complex balancing condition" for the general 
MAL [[621 [301. 

To formulate this condition for the MAL kinetics, let us start from the function G (12.141) and 
look for conditions that guarantee the inequality G < 0. 

For a given c°, we can rewrite the MAL reaction rate in the form 



( c \ a " 

(p r exp(a r ,V c G) = (prYly-^J , (2.16) 



Mv, (Pr ~ or)) exp[A(a r , V C G) + (1 - A) (ft, V C G)] . (2.18) 



where tp r = k r Yli(c^) a " = w r (c°). 

It is convenient to express G using an auxiliary function 9 of an auxiliary variable A: for any 
concentration vector, 

9(X) = Y,Vr exp[A(a r , V C G) + (1 - A) (ft, V C GQ] . (2.17) 

r 

This function is convenient because 

dfl(A) 
dA " 

r 

In this notation, 

G= -9\l). 

The function 0(A) is a sum of exponents. It is convex (9"(X) > 0). Therefore, if 9(0) = 9(1) 
then 9'(1) > 0. 

This condition, 9(0) = 9(1) (for all positive c) is called the complex balancing condition and it 
is sufficient for the dissipation inequality: 

G = -9'(1) < 0. 

The principle of detailed balance for the MAL equations has a form of existence of a spe- 
cial equilibrium point, a point of detailed balance. This existence implies important dynamical 
properties in the whole of Wt because of the very "rigid" monomial structure of MAL. 

The complex balancing condition also can be formulated as existence of a special "point of 
complex balance". Let us reformulate it in this way. 

Some vectors a r , (3 r for a given reaction mechanism may coincide. Let us denote by {y 1 , . . . , y q } 
the set of all different vectors a r , (3 r . For each yi we define Rf = {r \ a r = yi}, K[ = {r \ (3 r = 
yi}. In this notation, 
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^) = E E^( c °) exp(^V c G), 

) ( (2-19) 

= E I E ^ c °) J exp(2/i, V C G) . 

For any finite set of (different) vectors {y±, . . . , y q } the correspondent functions exp(?/j, V C G) 
of c G Wl are linearly independent because the Hessian of G is strictly positive definite. Therefore, 
the condition 0(0) = #(1) (for all positive c) is equivalent to 

£ w r (c°) = £ uv(c°) . (2.20) 

The point c° that satisfies (12.201 ) is called the point of complex balance and the complex balancing 
condition means that there exists such a strictly positive point of complex balance. In this case, the 
dissipation inequality, G < 0, with G defined by (12.141) holds for all positive points. 

Of course, a point of detailed balance is a point of complex balance as well. The reverse 
statement is not valid: for example, all the positive steady states of linear MAL kinetics (12.121 ) are 
the points of complex balance but they are not necessarily the points of detailed balance. 

The microscopic background for detailed balance is microreversibility, i.e. invariance of the 
microscopic classical or quantum equations with respect to the time inversion. 

The microscopic backgrounds for complex balance were formulated by Stueckelberg as unitar- 
ity of the S'-matrix II971 . It is necessary to add that the validity of the scattering (or Markov) model 
of elementary reactions is also needed: see [[541 and discussion in Section 1X1 



2.2. Mass Action Cell- Jump Formalism 

2.2.1. Stoichiometry of Diffusion Jumps 

We represent the physical space as a network of compartments. Each compartment is modeled as a 
cubic cell with an edge size I. The stoichiometric equations of diffusion describe interaction of two 
neighboring cells. To distinguish the quantities related to these two cells we use the upper indexes 
I and II (Fig. 0). 

The general stoichiometric equation for an elementary event of diffusion is (11.371 ) 

E <M + E E M + E ^ ■ ( 2 - 21 ) 

i i i i 

Coefficients a*'/ 1 , filf 1 are nonnegative. Usually, we assume that they are integers but in some 
situations real numbers are needed. 
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Elementary events (12.211 ) describe diffusion and do not include the transformation of compo- 
nents (reactions). Therefore, the total amounts of each component Aj coincide in the left and the 
right hand sides of ( 12.211 ): 

a l ri + a l H = + 0% . (2.22) 

Each elementary process (12.211 ) has two loss vectors, a]; 11 with coordinates a^j 1 and two output 
vectors, f3j ,u with coordinates /3^ n . Because of the conservation of particles of all types (12.221 ), the 
stoichiometric vectors of processes for the cells differ just by the sign of coordinates: 

Tr Tr Tr $r ' 

We define here a mechanism of diffusion as a system of stoichiometric equations for elementary 
events. The simple and basic examples are: 

• Fick's diffusion, A\ Af, Af -»■ Af, 

• The exchange of particles, A\ + Af ->■ Af + Af, 

• Clustering (diffusion with attraction), A] + sAf -)• (s + l)Af, sA\ + Af (s + 1)A\, 
s > 1; 

• Diffusion with repulsion, (s + l)A\ ->■ sA\ + Af, (s + l)A\ -)• A- + sA- 1 , (s > 0). 

Formally, diffusion with repulsion is the time-inverted process of diffusion with attraction (the 
porous medium model) but for s = 1 diffusion with attraction has no sense (the exactly uniform 
state cannot produce the nonuniform distribution). Therefore, the restrictions on s are different. 



2.2.2. MAL Equations for Diffusion 

Let us consider the system of stoichiometric equations (12.211) as a reaction mechanism for MAL 
(12.11) . If we apply MAL then the rate of the elementary diffusion process is 

W r ( C \ C U ) = k r 1[(C])< Y[(cf) a " . (2.23) 

i i 

For example, for Fick's diffusion, we have two elementary processes, A\ — > Af and Af — > Af 
The corresponding reaction rates are fcicj and k 2 cf. 

The composition of each cell is vector N I,IL . Components of this vector, Nj ,u = V 1 ' 11 ^' 11 are 
amounts of A; in the corresponding cell and V^ 1 ' 11 are volumes of cells. We describe the dynamics 
of the compositions of two cells by the equations: 

f4 = S ^rf/). (2.24) 

r 

where S 1,u is the area of the boundary between cells I and II. If there are many cells then 

(2-25) 
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with summation through all interacting pairs (I, J). 

For example, for Fick's diffusion, we have two elementary processes, A\ — > A] 1 and Af — > A\. 
The corresponding reaction rates are kic\ and k 2 cf . Equations ( 12.241 ) give 

= -S^hcl + S^k 2 cf . 

dt 

For several pairs, let us mention the symmetry in pairs: k\ = k 2 = k (we discuss this symmetry 
in more detail in the next subsection): 

For example, on a straight line (two neighbors), this equation gives 

^ = ks^(4 +1 + 4~ 1 -^). 

From this expression, the proper scaling of k with the cell size is obvious: 

dt V I 2 

the fraction c- +1 + c- -1 — 2c-// 2 approximates the second derivative and, hence, kS l ' 3 l 2 /V = 
const. 

For a cubic cell, V = SI and kl = const. 
2.2.3. Space Symmetry and Time Symmetry 

The system of elementary events should be symmetric with respect to space-inversion. For each 
elementary process (12.211 ) a space-inverted process is defined just by changing I to II and vice 
versa. We mark the quantities for the space-inverted processes by For example, 

i = -7 • 

Space inversion is an involution: if we apply it two times then we return to the original process. 

The key condition is: the rate functions for the space-inverted processes should differ just by 
transposition by vectors of variables, c 1 , c 11 (11.391 ): 

w' r {J, c 11 ) = w r (c u , c 1 ) . (2.26) 

For MAL this means that k r = k' r . 

The requirement of space symmetry distinguishes diffusion from various types of advection and 
transport driven by external force. This condition is necessary for existence of diffusion equations 
when the cell size / — > (see the next subsection). 
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Inversion in time differs, in general, from the inversion in space. For example, for the ele- 
mentary process 3A 1 — > 2A l + A u the space inversion gives 2A n — > A 1 + A u (we exchange the 
upper indexes, I— >ll and II— >T), and the inversion of elementary events (T-transformation) gives 
2A 1 + A 11 3A 1 (here, we change direction of arrow). We have to stress that inversion in time 
assumes micro-reversion. At the macroscopic level it does not mean change t to —t in the kinetic 
equations but the transformation of the direct processes into reverse ones (inversion of collisions, 
for example). 

Time symmetry (microreversibility) means that the principle of detailed balance is valid. In 
this case, all the consequences of the principle of detailed balance are applicable (Section \2. 1.3.1) . 
a global Lyapunov functional exists, and every positive equilibrium is a point of detailed balance. 

Microreversibility and symmetry of space are independent properties of the system. Never- 
theless, for some elementary processes, the space-inverted process coincides with the reverse (the 
time-inverse) process. If a diffusion mechanism is constructed from such processes then the sym- 
metry in space is equivalent to symmetry in time (to the principle of detailed balance). 

This specific class of diffusion mechanisms includes such mechanisms as Fick's diffusion or 
the diffusion by exchange of particle positions. 

Indeed, for Fick's diffusion, when we exchange the upper indexes in an elementary process 
A 1 — > A 11 (inversion in space) then we get the reverse process as well, A u — > A 1 (inversion of 
arrows). Analogously, for A 1 + B n — > A 11 + B 1 inversion of space gives the reverse process as 
well: A u + B l -> A 1 + B n . 

This fundamental property is formulated in the following theorem. 

Theorem 2. Let a complex diffusion process consist of elementary processes, which satisfy the fol- 
lowing property: the space-inverted elementary process coincides with the inverse process. Then, 
for the mass action law equation of diffusion A2.25\) . the principle of detailed balance is valid, the 
global convex Lyapunov functional exists and the uniform distribution is asymptotically stable. 

Indeed, due to space symmetry, a uniform distribution is an equilibrium and each process is 
equilibrated at this state by its space-inverted process. At the same time, this distribution is a point 
of detailed balance. Therefore, the results about the principle of detailed balance are applicable. 

2.2.4. Arrested Diffusion and Boundary Equilibria 

Existence of a uniform distribution, which is a point of detailed balance, existence of the global 
Lyapunov function and asymptotic stability of the uniform equilibrium distributions do not mean 
that there exist no nonuniform equilibria. 

The phenomenon of boundary equilibria is well known. For example, an autocatalytic re- 
versible reaction A + B ^ 2 A has two equilibria for a given value of a stoichiometric conservation 
law, b = ca + cb =const. One equilibrium is strictly positive, with positive concentrations of A and 
B. Another is a boundary equilibrium, ca = 0, eg = b. The positive equilibrium is asymptotically 
stable, the boundary equilibrium is unstable but if the initial state is near the boundary (ca is close 
to zero) then slow relaxation occurs [|43l . and the motion may be arrested for a long time near this 
state. 
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There are well known effects of arrested diffusion caused by changing of temperature. The 
solid solutions show the effects of diffusion, which has been arrested by chilling below a threshold 
temperature in a very short time [|74ll . 

Kinetic effects of arrested diffusion are also possible. For example, let us consider the diffu- 
sion mechanism by jumps to free places (Fig. |3(b)[ ): any distribution of components Ai for zero 
concentrations of free places Z is stationary, and for a small concentration of free places diffusion 
is slow. 

For effects of arrested diffusion, the average concentration of free places should not be small. 
There may be, for example, a layer of particles with low mobility between a dense island of parti- 
cles with high mobility and an island of free places. Formally, it is possible to construct many such 
situations. All of them may be characterized as follows: either a small change of concentrations or 
a small change of constants (or both) lead the system to a nonuniform equilibrium. At this nonuni- 
form equilibrium some of the concentrations in some cells take zero values and, therefore, some of 
fluxes are also zero. Because of the appearance of these zero concentrations, such an equilibrium 
is called a boundary equilibrium. 



2.3. Continuous Diffusion Equation 

2.3.1. MAL Diffusion Flux 

Let us consider an elementary process together with its space-inverted process 

E a «4 + E ^ E fcA + E Prl^Y , 

i i it 

E <44 ? + E «M E ^ + E fcti ■ 

i i i i 

The reaction rates are 

Wr(c\c n )=k r ]J(cl)<]J(c?)<, 

i i 

w' r {c\ c n ) = w r (c l \ c 1 ) = k r H(cf)< H(4) a " , 

i i 

where we take k' r = k r due to the symmetry in space. 

To first order in /, the flux vector for Ai in this process is 

Jri = -lri[w r {c{x), c(x + /)) - W r (c(x + I), c(x))] 

dw^c 1 , c 11 ) dw r (c l ,c 11 ) 



dc 1 } 



-l^ ri w r (c(x),c(x)) E 



C I =C II =C (^ 

a 11 - a 1 ■ 
— -Vc,-(x) 



del 



V Cj (x) 



c I =c II =c(x) j 



C 



-ik lri n^ 9+< e 



(2.27) 



(2.28) 



(2.29) 
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Here, 7 ri = (5 l ri — a l ri (input minus output in the first cell); the minus in front of the formula appears 
because the direction of flux from cell I to cell II (from x to (x + I)) is positive. 

The factor 1/c, never leads to a singularity in the flux because cj enters in the monomial 

Ylq Cq rci+arq with the power a\j + a". This power is strictly positive if the coefficient (a" — a^) 
is not zero. 

The proper scaling of k for grid refinement or coarsening is kl = d = const in order not to 
change the first order expression for flux (I2.29I ). 

According to (12.291 ), the matrix of diffusion coefficients for the elementary process (12.271) (to- 
gether with its space-inverted process) is 

D rij ( C ) = d [n c « r * j ^ — ' (23o) 

where d = const (= kl). 

The corresponding diffusion equations have the divergent form: 

3c 

— = div(D(c)Vc) , (2.31) 

where c is the vector of concentrations and D is the matrix of diffusion coefficients (12.301) . 

It might be useful to represent the flux ( 12.291 ) similarly to the Teorell formula (11.71) . For this 
purpose, let us collect under V the terms which represent the chemical potential in perfect media: 
/j, = RT In c + fj, . We assume that T and fx are constant in space. With these conditions, 

Jn = ( II ^ q+< ) £(°S - VA*i(x) • (2.32) 



2.3.2. Examples 



Let us illustrate application of formula (12.291) by several elementary examples. 

First of all, Fick's law: A 1 -> A u and A u -> A 1 . For this system, a 1 = 1, a 11 = 0, /3 1 = 0, 
P n = 1 and 7 = -1. 

Formula (12.291 ) gives 

J = Ikc—Vc = -IkVc. 

c 

This is exactly the standard Fick law. The diffusion equation is d t c = dAc. Here and further in 
this subsection we use d for Ik. 

Exchange of positions: A 1 + B u — > A 11 + B l together with the space-inverted process A 11 + 
B l — > A 1 + B u , which is the same as the reverse process. For this case, a\ = l,a l B = 0, a l \ = 0, 
1, p A = 0, p\ = 1, 0° = 1, 0% = 0, lA = p A - a\ = -1, and lB = f3 l B - a l B = 1. Due to 



a 



B 



Ja 
Jb 



-d(-l)c A c B 
-dc A c B 



— Vc A + — Vc B 

c A c B 



-1 1 „ 

— Vc A + — Vcg 

ca c b 



-d(c B Vc A - caVc b ) 
d(c B Vc A - c A Vc B ) . 



(2.33) 
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The diffusion equations are 

d t c A = d(c B Ac A - c A Ac B ) , d t c B = d(c A Ac B - c B Ac A ) ■ 

Repulsion of components A and B {A is mobile). The mechanism is A 1 + B l A 11 + B l 
together with the space-inverted process A 11 + B u — > A 1 + B u , which does not coincide with the 
reverse process. For this mechanism, a l A — 1, a B — 1, a 1 } = 0, a" = 0, (5\ = 0, /3 B = 1, (5\ l = 1, 
(3 l B = 0, 7,4 = f3 A — a l A = —1, and 7b = (5 l B — a l B = 0. Formula (12.291) gives: 



J A = -d(-l)c A c B 



—Vc A + — Vc B 

c A c B 

(2 34) 

-d{c B Vc A + c A Vc B ) = -dV{c A c B ) , 



Jb = 0. 

We can see that B activates diffusion of A (the term c B Vc A ) and, at the same time, pushes A in 
the area with lower concentration of B (the term c A Vc B ). 
The diffusion equation is 

d t c A = dA(c A c B ) . 

The no-flux steady states of these diffusion equations are given by the condition: c A c B = const. 
If we assume that B is also mobile by a similar mechanism then we get a system of equations 
(with two different diffusion coefficients): 

Ja = -d A V{c A c B ) , J B = -d B V{c A c B ) , (2.35) 
d t c A = d A A(c A c B ) , d t c B = d B A(c A c B ) . (2.36) 



Let us change variables: 



d A d B d A d B 



In these variables, c A c B = (c 2 + — c 2 _), <9 t c_ = and for c + we have the porous media equation: 

d t c + = ^Acl . 

2.4. Principle of Detailed Balance and Dissipation Inequality 

2.4.1. Detailed balance and Coupling of Direct and Reverse Processes 

In this subsection, we formulate the principle of detailed balance for MAL diffusion. Physically, it 
follows from microreversibility. 
For every elementary process, 

]T a\ t A\ + £ c^Af ]T P l ri A\ + P'AY , (2-37) 

i i i i 
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the reverse process is (just invert arrows): 

E pia + ASM? -+ E + E a " A " • < 2 - 38 ) 

i i i i 

Let us distinguish the quantities for the reverse and direct processes by the upper indices ±. The 
simple algebraic relations hold: 



1,11+ fll.ll* d 1,11+ 

Therefore, 



off + = #T + = °t + «S" = & + + ^ + • ( 2 - 39 ) 
The reaction rates are 



(2.40) 



Let there exist a uniform strictly positive point of detailed balance: such a strictly positive 
vector c* that for all r 

Wr[C ,C ) = W r {C ,C ) . 

This means that 

i i 

Let us use the relations a l ri + a" = + Therefore, the principle of detailed balance for 
MAL diffusion can be reformulated: for all r 

k r — — k T • 

Let us join processes (12.371) . (12.381) and write for them the diffusion flux analogously to ( 12.291 ): 

Jri = -Jri[Wr(c(x), C(X + I)) - W~ (c(x + I), c(x))] 



-hri E 



dw+(c x ,c 11 ) 



dcf 



dw r (c 1 , c 11 ) 



dc l 

c i =c ii =c(a .) ^ 



c i= c ii= c ( x ) / 

IrlxirJ^ (2.41) 



= -/7 ri u7 r (c(x), c(x)) ^ — Vc^x) 

Co 

\ q / j i 

According to this expression for the diffusion flux, (12.411) . the matrix of diffusion coefficients 
for the elementary process (12.271) (together with its space-inverted process) is 

Dljic) = d r foj c < +a A l^Ti , (2.42) 
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where d = const (= k r l). 

This matrix is symmetric with respect to the inner product 



b) c = J2 — ■ (2-43) 



This means that 

(Da,b) c = (a,Db) c 

or in coordinates 

Dijdjbi _ a 



ij 



Indeed, let w r = d r (j^^ c ^ rq+arq \ _ p or two scalar products, (Da, b) c and (a, Db) c we get 

(Da, b) c = = «V E 7 " 7n ' aA 



and 



/ ST^ CliDijbj _ sr-^ ai^ r i^ r jbj 

(a, Db) c = ^ — = ^ / 



These two expressions differ only in the notation of dummy indexes, i and j. Therefore, D is 
symmetric with respect to the inner product (12.431) . 

D is also positive semi-definite. Indeed, for any vector £, 

<^O = £^ = ^,7) 2 >0>. 

i Ci 

This expression may be zero at a positive state (q > 0) if and only if the vector £ is orthogonal 
to the vector 7 r in the inner product (12.431) . If we summarize diffusion coefficients for all pairs of 
mutually inverse elementary processes then we get 

^ = Ew = E^- : t m - (2 - 44) 

r r j 

For this D^-, 

Q 

and is zero if and only if vector £ is orthogonal to all vectors "f r in the inner product (12.431) . 
The corresponding diffusion equations have the divergent form: 

Be 

— = dw(D(c)Vc) . (2.45) 
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It might be useful to represent the flux (12.411 ) similarly to the Teorell formula ( 11.71 ). For this 
purpose, let us collect under V the terms which represent the chemical potential in perfect media: 
fx = RT In c + fx . We assume that T and fx are constant in space. With these conditions, 

Jri = -^7n ( J] °f q+< J E 7riV^(x) . (2.46) 

The difference from the MAL Teorell formula (12.321) is obvious: for the systems with detailed 
balance (12.461) the matrix of the coefficients in the Teorell formula is symmetric for each elementary 
process together with its reverse process because it is a product of a number (in square brackets) 
and the symmetric matrix ^ r Cirr 

1 U I T TT \ 

Irilrj ■ (2.47) 

For the general MAL diffusion (with detailed symmetry in space) the matrix these coefficients 
for the elementary process together with its space-inverted process is not symmetric in general: it 
is a scalar multiple of the matrix 

Its symmetry may be guaranteed if the space inversion of the elementary process coincides with 
its reversion in time (i.e. a* 1 coincides with /?*). 

For MAL chemical kinetics, there is a sufficient algebraic condition for detailed balance that is 
independent of the microreversibility and follows just from the stoichiometric equations. Indeed, 
let us assume that all the reactions are reversible and all the stoichiometric vectors 7 r in (12.71) are 
linearly independent. Then, at the equilibrium, from the condition c = ^ r w r j r = we get 
w r = for all r. That is the detailed balance condition. 

For MAL diffusion we have also a specific ("diffusion") algebraic condition: if for all elemen- 
tary processes the space-inverted reaction is the reverse reaction then the principle of detailed bal- 
ance follows from the space symmetry condition. For such systems, the coupling "process-space 
inverted process" coincides with the coupling "process-reverse process" and equations (12.411) co- 
incide with (12.291) . 

For an elementary process (12.371) . let the space-inverted process not coincide with the reverse 
process. With the condition of detailed balance, it is straightforward to check that for the pair of 
elementary processes (12.371 ), ( 12.381) , the space-inverted process to (12.371) and the reverse process to 
this space-inverted one generate the same diffusion equations and the diffusion coefficient as the 
original pair does. These two couples together produce the flux that is twice as large as (12.4- lb . 

2.4.2. The Dissipation Inequality and Detailed Balance 

In this subsection, we construct the Lyapunov functional for the general MAL diffusion models 
and calculate its time derivative due to the diffusion equations. 
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Let us select any strictly positive reference concentration vector c* and take (12.8 



c 



Let us consider this system in a bounded domain K with smooth boundary and with zero fluxes 
through its boundary: (n, J) = at any point of dV at any time (n is the vector of the outer 
normal). Due to the definition of flux (12.291) . it is sufficient that (n, grade.,) = on dV for all j 
(but it is not necessary). 

The Lyapunov functional is 



G = £ 



dx . (2.48) 



Due to the boundary conditions and the Gauss-Ostrogradskii theorem 

f = - £ X k ( |) ^ = E / v (V. (>* |) ■ *) ■ (2.49) 

Let us assume the principle of detailed balance and calculate G (12.491) due to diffusion equation 
with the matrix of diffusion coefficients (12.441) : 



(2.50) 

dx 



= _ 5Z/ ^r^(V x (lnc i ),7 r i7 n V a; (lnc i )) dx 
= -^2 J w r ^V x ^2 ln (7«Cj) j , V z In (7 rj -c 3 - 
= ~J2 [ ^r(V x .( 7 „V c G)) 2 dx<0. 

r 

Here, ( , ) is the standard Euclidean inner product both in space and in the concentration space. 

As we can see from this dissipation inequality, the only positive equilibria for diffusion equa- 
tions with detailed balance conditions satisfy the conditions (7 r , V C G) = const for all r. Another 
form of these conditions is 

J cj" = const . 

i 

Inequality (12.501) demonstrates a significant difference between the two classes of diffusion 
mechanism. If the stoichiometric vectors 7 r form a basis in the concentration space then all the 
equilibria are uniform because in this case the condition (7 r , V C G) = const (for all r) implies 
V C G = const, that is, In q = const for all i, hence, q = const for all i. Let us call this mechanism 
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the mechanisms with mixing. The first example is Fick's mechanism: if the diffusion constant is 
not zero for all components then all the equilibria of the system are uniform. 

If span{7 r } does not coincide with the concentration space then there exist the invariants of 
diffusion. They are given by the linear functionals that annihilate all the vectors j r . For example, 
in diffusion by jumps on the free places (Fig. |3(b)[ ) the value of the sum cz(x) + J2 i Ci{x) is locally 
conserved. In the mechanism, A 1 + B l — > A 11 + B 1 , A 11 + B 11 — > A 1 + B n , concentration of B 
is conserved. These locally conserved quantities together with the condition of positivity q > 
define a convex body where the vector of concentrations may be situated at a given point. This 
body depends on the values of the conserved quantities, differs for different points, but does not 
change in time. 

2.5. Complex Balance in MAL Diffusion 

2.5.1. Complex Balance Conditions for MAL Diffusion 

The complex balance condition does not assume any space or time symmetry. The only micro- 
scopic assumption is the Markov fast microscopic kinetic with relatively small amount of active 
intermediate complexes ||97ll54l . 

We discussed this condition for the MAL kinetics in Section 12. 1 .4.1 and now let us transform 
it into a condition for the MAL diffusion equation. First of all, we should abandon the symmetry 
conditions k' = k (space symmetry) and k + = k~ (microreversibility). Without these conditions, 
the zero-order terms in the expression for fluxes will not be annihilated by balance between the 
elementary processes with given pair of vectors (a£, a* 1 ) and the processes with the same pair 



where all the elementary processes have different numbers r and the space-inverted and reverse 
processes are represented separately. Let us consider all pairs of vectors, (a*, a* 1 ) and 




E + E E + E 



ii 

i i 



(2.51) 




The complex balance condition is: there exists a strictly positive vector c* such that 




(2.52) 



reRt 



r£R q 



For MAL this means 



E k r\\{c*T l+aV = E k i\{«T lq+al " for all/. 



(2.53) 
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For r G R q , 



al + a? = 2/J + y? , 



and for r G R t 

Pl + f3l l = yl + yf. 

Therefore, all the monomials coincide in the right and the left hand sides of (12.531) for the given q 
and we can write the complex balance condition for diffusion in the following form: 

k r = K for all q . (2.54) 

Let us calculate the flux 

J = 7 r uv(c(x), c(x + /)) 

r 

at first order in I. We group terms in this sum and each group will corresponds to a pair y q . For each 
elementary process with number r there are two q, q a (r) and qp{r): r G R qa ^ and r G R q ^ r y We 
split the term ^ r w r (c(x), c(x + /)) in two terms: 

j r w r (c(x), c(x + I)) = (3 r w r (c(x), c(x + I)) — a r w r (c(x), c(x + /)) 

We associate the first of them with y qf) ( r ) and the second with y qa ( r )- The result is represented 
below: 

J = — Y^ 7rW r {c(x), c(x + /)) 

r 

= 2^ a l w r(c(x), c(x + I)) — ^ 0lw r (c(x), c(x + I)) 

1 \r€R+ r<=R- 

= y l\ W r{ c { x ), c ( x + 0) - W r-(c(x), c(z + I)) 

1 \rei4 r€ 

+ ' E i4 II ^ E^T^-E^-^h-c) 



rdRn 



(2.55) 



C,' 14 — ' 14 — ' c 



The zero-order term is zero because of the complex balance condition (12.541) . Let us take into 



account that both for r G R q and r e R q 



a 1 + a 11 = y\ + y\ l 
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Therefore, if for two y q , y s 
then 

y\ + yf = y\ + yf ■ 

For the set of all vectors y l q + y q l we use notation Z. For each z E Z 

Y z = {y q | y\ + y\ l = z} 

Let us group the terms with the same vectors z = y l q -\-y l q E Z together and write the expression 
for the flux for the first order: 



a?, ^ o& 



(2.56) 



Vi&z 3 \r£R+ 3 r&R- 3 J 

Let us study the expression for J z for given z. First of all, 

Vq&z j \r€R+ 3 r£R~ 3 j 

because each k r from this sum enters it twice, with the opposite signs, one time as an element of 
Rt for some y q E Y z (with +) and the second time (with — ) as an element of R~ for another 
y s E Y z with the same sum 

y\ + yf = y\ + y l q l = z. 

We note that y l = z — y l q for y q E Y z in the sum ( 12.561) . This allows us to rewrite the expression 
for J z : 

= " E »? E ( E *vf " E ^ J Vc, • (2.57) 

ViZYz 3 \r£R+ 3 r&R- 3 J 



2.5.2. The Dissipation Inequality and Complex Balance for MAL Diffusion 

We would like to demonstrate some similarity of the expression for J z and some formulas from the 
theory of Markov chains. Vectors y q E Y z numerate the states. Elementary processes correspond 
to transitions between states. Each nonzero constant k r corresponds to two vectors y q , y s E Y z : 
r E Rt and r E R~. We substitute the index r by two indexes q, s and use notation k sq (or even 
k s ^. q ). If there is no nonzero constant for this pair q, s then we take k s <_ q = 0. In particular, 
k qq = 0. The complex balance condition (12.541) reads: 

k ^ = J2 k os- (2-58) 

9 1 



47 



A.N. Gorban et al. 



Quasichemical Models of Multicomponent Diffusion 



This means that the constants k sq describe a continuous time Markov chain with the Master 
equation 

Tr q = y^(fc gs 7T s - k sq n q ) (2.59) 

s 

and equidistribution in equilibrium. Here n q is the probability to find the system in the state q and 
k qs ir s is the probability flux from the state s to the state q. We can use the steady state condition 
(12.581) and rewrite the Master equation (12.591 ): 

Kq = £ Ksi^s ~ TTg) • (2.60) 

s 

From this form, it is easy to see that the functional 



2 ^ q 



q 

monotonically decreases due to the system dynamics: 

— = 22 ^g k qs(^s - 7T g ) < . (2.61) 
sq 

To prove (12.611) . let us use the identity 

£ 7Tg = y^(fe gs 7T 6 - k sq TX q ) = ^ Ksi^s - 7T g ) = . (2.62) 

q s s 

This condition holds for all values of numbers n q (this is obvious for the Master equation in the 
form (12391) ). 
In particular, 

s ' 

Let us add this expression to the right hand side of (12.611) : 

because > and the expression in the parentheses is non-positive: 

11 1 

Therefore, we have proved the inequality for our set of coefficients k qs : 

£ n q (k qa n a ~ k sq -K q ) = ^ Kqkgsi-Ks - n q ) < (2.64) 

sq sq 
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for any set of numbers ir q . This inequality means that dH/dt < (12.6 1 1 ). It is zero if all ir q coincide 
(tt.s = 7r g for all s, q). 

For our purposes, it is important to know when the zero time derivative of H (dH/ dt = 0) is 
equivalent to the equidistribution (n s = n q for all s, q). They are equivalent if the Markov chain 
(12.591) is ergodic. The conditions of ergodicity are well known 11951110611 : the chain (12.591 ) is ergodic 
if for any two s, q (s ^ q) there exists a oriented path from s to q in the graph of the network, that 
is such a sequence 



This means that the graph of of transitions of the Markov chain (12.591 ) is strongly connected. 

To apply this inequality to the proof of the dissipation inequality, we have to rewrite the ex- 
pression for J z (12.571 ) using these notations, k qs instead of k r : 



Now we are in the position to prove the dissipation inequality for MAL diffusion equation 
with complex balance. In a bounded domain V with smooth boundary and without fluxes through 
boundary we have to estimate G (12.481) . (12.491) . 



7*0, ri, . . . , r 9 that s = r , q = r g and k, 



> Ofoallj =0,... 




(2.65) 



Vq& z 



S 




(2.66) 



Due to the representation of J z (12.651) . 



^(V.lnc,-, J/) = J2^( k ' 



•qs 



n s - k sq it q )) < , 



(2.67) 




where n q is a space vector: n q = J2i y l qiS x In Q and for y q , y s 

vis + y l q, s = z ■ 



This expression has exactly the form (12.641) for each space coordinate. Finally, 




(2.68) 



and it is zero if all 7i q = y t 



\V X In Q coincide. 
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The reverse statement, 

all Ti q = ^ In Cj coincide if ^(V x In cj, Jf) — , 

* j 

is true if the auxiliary Markov chain is ergodic for given z (i.e. the graph of transitions is strongly 
connected). Let us assume this ergodicity. For every z we can define a linear subspace E z in the 
concentration space given by the system of equation 

E z = {e\ (j/J, e) coincide for all y q e Y z } . 

If 

n & = {o} 

zez 

then all the equilibria for this mechanism of diffusion are uniform. In particular, they are uniform 
if at least one y q \ = 0. 

2.6. Intermediate Summary 

We presented the formal language of the stoichiometric mechanism for description of nonlinear 
diffusion (12.211) . (12.271) . (12.371) . (12.381) . The general construction of the diffusion equations under 
various conditions is given: 

• For systems with symmetry with respect to inversion in space (k r = k' r ) (12.291 ). (12.301) ; 

• For systems with microreversibility (k+ = k r —) (12.411) . (12.421) ; 

• For general systems with Markov microkinetics, which satisfy the complex balance condi- 
tions (12341) . (12361) . 

For systems with detailed balance (microreversibility) and with complex balance (Markov mi- 
crokinetics) the explicit formula for the free energy-type Lyapunov functional G (12.481 ) is 



«M3 - 1 + E« 



dx . 



We found that G < for the systems with detailed or complex balance (12.491 ). (12.661) . These 
inequalities guarantee the thermodynamic behavior of diffusion. The detailed symmetry with re- 
spect to the inversion in space is insufficient for such an inequality and the diffusion collapse for 
them is possible. 

The Mass Action Law by itself does not imply thermodynamics. In that sense, it is too flex- 
ible and needs additional requirements to respect the basic physics. This redundancy of MAL 
allows, at the same time, to use them in many other areas. The famous Lotka and Volterra models 
in Mathematical ecology J73l 1 1 10H are implementations of MAL for description of surviving and 
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reproduction of animals, that is far from the initial area of MAL applications. Creation of mathe- 
matical genetics [33 ] and analysis of dynamical aspects of biological evolution [|38l also use MAL 
in their backgrounds. Combination of MAL kinetics with nonlinear diffusion is very important 
for analysis of biological invasions and other phenomena in ecology 11571 |6T1 [%1 [89]|. The MAL 
for diffusion can also generate equations that have no direct physical sense but may be used for 
modeling of some phenomena of non-physical nature. 

On the other hand, in this approach, we did not take into account some basic physical proper- 
ties, namely, the momentum and the center of mass conservation. The diffusion transport should 
be coupled with the viscous transport or elastic deformation (or both) because two reasons: 

• The mass average velocity of diffusion 

Y,i m i J i 

Li m i 

where m t is the molecular mass of the A { particles, is, in general, not zero; 

• The change of the mixture composition implies the change of pressure and, hence, the vis- 
cous flux or the elastic deformation (or both). 

The careful analysis of these effects should give, for example a theory of the Kirkendall effect. In 
1942, Kirkendall demonstrated experimentally that different atoms can migrate at different rates in 
an alloy, and this diffusion is accompanied by measurable local volume change and displacement 
of interfaces |j67ll8TI . The kinetic theory of this effect is still under development and there remain 
open problems and new ideas are needed ll82l . 

We return to the problem of correct description of the general transport equation in next Section. 

The MAL formalism for diffusion is a flexible and effective tool for modeling. The semi- 
discrete MAL model may be used for numerical modeling directly as a sort of finite elements. 
Their coarse-graining and refinement should follow the main rule kl = d where I is the cell size, 
A; is a kinetic constant for the finite elements, d is the invariant diffusion coefficient. For unstable 
processes, these provide a biharmonic regularization. These kinetic finite elements respect the 
basic physical properties like positivity of concentration, conservation laws and the second law of 
thermodynamics (under the relevant conditions of detailed or complex balance). 

3. Generalized Mass Action Law for Diffusion 

3.1. Free Energy, Free Entropy, Chemical Potentials, 
Activities, and Generalized Mass Action Law 

3.1.1. Thermodynamic Potentials 

In this Subsection, we present the thermodynamic approach to the generalized MAL. Exactly as 
in Section 12.1.1 we start from the chemical kinetic equations and then extend our approach to the 
transport processes. 
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In the thermodynamic approach, the kinetic description of the multicomponent system requires 
the following inputs 

1 . A list of components; 

2. A thermodynamic potential; 

3. A list of elementary reactions represented by their stoichiometric equations; 

4. A set of reaction rate constants. 

Exactly as it was in Section l2TTl the list of components is just a set of symbols (the component 
names). We usually assume that this set is finite, A t , A 2 , . . . , A n . The definitions of stoichiometric 
equation and the corresponding vectors a ri , f3 ri are also the same. 

There are many thermodynamic potentials and they form two series: energy and free energies 
and, on another hand, entropy and free entropies (the Massieu-Planck functions). Each of them 
has its own "natural variables" and if one of them is given in the natural variables then all other 
thermodynamic functions can be produced. 

The names and standard notations of variables are: 

• Internal energy, U ; 

• Entropy, S; 

• Enthalpy, H; 

• Temperature, T; 

• Volume, V; 

• Pressure, P; 

• Number of particles (or moles) composing the ith component A4, iVj (N is vector with coor- 
dinates Ni, the vector of composition); 

• Chemical potential of the zth component A±, fa. 

The first potential in the energetic series is the internal energy U (S, V, N) in the natural coor- 
dinates S, V, N and 




The enthalpy, H(S, P,N) = U + PV has the natural coordinates S, P, N and 



dH = TdS + VdP + Vi dN i ■ 
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The free energy (the Helmholtz energy), F(T, V, N) = U -TS and 

dF = -SdT - PdV + VidNi . 

i 

The free enthalpy (the Gibbs energy), G(T, P,N) = H -TS and 

dG = -SdT + VdP + Vi dN i ■ 

i 

The grand potential, Sl(T, V, //) = U - TS - £\ f^Ni and 

dVl = -SdT - PdV - N ^ ■ 

i 

The entropic series starts from the entropy S(U, V, N) and 

dS = LdU+^dV-Y §diV t . 

rp rp / j rp I 

i 

Therefore, the main set of the intensive variables for the entropic series is 

1 _ dS(U, V, N) P _dS(U 1 V,N) _ dS(U,V,N) 
f ~ dU ' T ~ dV ' ~f~ ~ dNi ' 

The Massieu function, $(T _1 , V, N) = S - T~ X U (= —F/T) and 

d$ = -ud (I) + ^dv-J2 j m . 

The Planck function, -(T^ 1 , T^P, N) = S - T~ X U - T^PV (= —G/T) and 

d* = -w(i)-vdg)-j:fdJVi. 

All these functions are used for the definition of equilibrium. The main definition for an isolated 
system follows the R. Clausius two main laws formulated in 1865 [TTTl : 

1 . The energy of the Universe is constant. 

2. The entropy of the Universe tends to a maximum. 

More precisely, the entropy of the isolated system tends to a maximum under given U, V and values 
of other conservation laws. Let the conservation laws be given: 

i 
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then the equilibrium is the maximizer of the entropy under given values of U, V and Bf 

S(U, V, N) -> max subject to given U, V and {Bj} . (3.1) 

Gibbs ll39l paid much attention to the dual formulation of this condition: 

U(S, V, N) ->■ min subject to given S, V and {Bj} . (3.2) 

Other definitions of the same equilibrium are available through the free energies and entropies. For 
the free energies this is condition of minimum. Analogously to (13.21) we get 



H(S, P, N) -»■ min subject to given S, P and {Bj} , 

F(T, V, N) min subject to given T, V and {Bj} , (3.3) 

G(T, P, N) -> min subject to given T, P and {Bj} . 

For the free entropies the equilibrium should be the maximizer: in addition to (13.11 ) 

( \ \ 1 
$ I — , V, N J — > max subject to given — , V and {Bj} , 

(3 4) 

/l P \ IP V ' 

S I - , — , JV J -»■ max subject to given -, — and {Sj} . 

In extensive thermodynamics, V, U, S, and are the extensive variables, that is they are 
directly proportional to the system volume if we just join several copies of the system with the 
proportional increase of the volume. Therefore, U (S, V, N) is the homogeneous function of the 
first order, and the equation for dU can be easily integrated (this is the Euler theorem): 



U = TS - PV + J2^ N i- 



In this case, 



H = U + P V = TS + 2_^ ViNi , 

i 

F = U _ TS = -PV + J2^ N i , 

i 

G = H-TS = J2^ N *> 

i 

n = u -ts -J2 = - pv > 

i 

Free energies have a very important physical and technical sense. They measure the available 
work under given conditions. 

Free entropies coincide (up to some constant additions) with the entropies of the minimal iso- 
lated system, which includes the system under consideration. This statement was analyzed in detail 
in 11421 but is still not very well known. Therefore, let us prove it. 
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Physically, when we consider a system under isothermal condition, this means that the system 
is in contact with a large thermal bath. The state of the thermal bath is characterized by two 
variables, Ut and Vp- The entropy of a thermal bath is St(Ut,Vt). The total entropy of the 
isolated system "a system + the thermal bath" is 

S(U, V, N, U T , V t ) = S(U, V N) + S T (U T , V T ) ■ 

The equilibrium is the maximizer of the total entropy S for given total energy U = U + Ut, 
values of volumes V, Vt and linear conservation laws {Bj}. In particular, 

d[S(U,VN) + ST(U-UVT)} _ Q 

dU ' (3 5) 

dS(U,V,N) dS T (U T ,V T ) 



dU dU 



u T =u-u 



This means that the temperatures of the thermal bath and the system are equal, T = Tt. 
Let us take the conditional maximum function (under conditions U + Ut — U, T = Tt'. 

S*(U, V N, V T ) = S(U, V N) + S T (U - U, V T ) 

= S(U,VN) + V T S T (^-^i 1 

(U \ U (3.6) 
= S(U,V,N) + V t St f — ,1 j - — 

= $ + VtSt +0(VC [ r 1 ) 

This function differs from the free entropy $ by a constant St(U, Vt) and an infinitesimal 0{V^ X ), 
which goes to zero when the bath increases. 

Therefore, the free entropy \& (the Massieu function) is equal to the entropy of the minimal 
isolated system, which includes the system under consideration and the large thermal bath (up to a 
large constant and an infinitesimal additions). 

For the Planck function S we have to consider a system under a constant pressure in the contact 
with the same large thermal bath. The only difference is that instead of the internal energy of our 
system we have to take the enthalpy. The enthalpy is the energy of the system plus the device, 
which keeps the pressure constant (potential energy of a heavy piston of the given weight). So, the 
total energy is the energy of the minimal isolated system U — H + Ut and everything else is the 
same as for $: the free entropy is the entropy of the minimal isolated system, which includes the 
system under consideration, under condition of the partial equilibrium with auxiliary systems and 
up to a constant summand. 
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For perfect systems, by definition, 



i 

PV = RT ^ Ni , 



(3.7) 



where Uj(T) is the energy of one mole of A; at the temperature T. 

Under this assumption, the entropy S is defined up to an arbitrary uniform function of first 
order S (N): 



-R In c, 



S = RS (N) + J2 N i 

i 

If we assume that S (N) is a linear function, 



T 



To 



r dr 



dr 



(3.8) 



then 



— R(\n Ci — 5i 



T 



1 dui(r) 
t dr 



dr 



(3.9) 



where q = N/V is the concentration, T > is a reference temperature (we assume that on the 
interval [T , T] the system is perfect (I3.7|) ). 

It is necessary to stress that linearity of S (N) does not follow from the assumption (13.71) and 
is an additional hypothesis. For a general perfect system (13.71) we can state that S (N) is a uniform 
function of the first order only. 

Formulas (13.71) . (13.91) allow us to express the free energy in the proper variables, F(T, V, N): 

F(T, V,N) = U-TS 



^N iUi {T) + BT^Ni 



In ( ^ I - 5, 



T 



1 duj(r) 
Rr dr 



dr 



(3.10) 



From this formula for F(T, V, N), all other thermodynamic functions can be expressed locally: 
OF n _ F _ T dF OF _dF_ 

dT ' &T ' d V ^ dNi ' " " " 

Generalizations of the free energy for non-perfect systems are often produced by transforma- 
tions of (13.101) . The first generalization describes a system of small admixtures to a general system. 
Let the "background" system have the extensive state variables M. Interaction of small admixtures 
is negligible and the formula PV = RT Ni describes the osmotic pressure of the admixtures.. 
Then we can write, analogously to (13.101) : 



F(T, V, N, M) =F (T, V,M) + J2 N * u * ( T ' y) 

+ rtJ2 n < 



Ni 
V 



To 



1 dUj (r, f ) 
Rt dr 



dr 



(3.11) 
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Here, the energies Ui and parameters Si are functions of densities M/V. For each given value of 
these densities, this formula coincides with (13.101) . 

The first model of non-perfect gases is the van der Waals gas. To write the free energy for this 
type of gas (or gas of admixtures) we have to take into account two effects: the excluded volume 
per mole of particles Ai,Vi, and the energy of attraction for particles Ai, Aj with the energy density 



(Xij C{ Cj 



(minus because this is attraction). 

For the free energy these effects give: 



F(T, V,N)=J2 N ^T) ^44 

i 

+ rtJ2 n > 



111 



v v 

Ni 



V-E t ViNi 



Si 



1 dui 



To 



Rt dr 



dr 



(3.12) 



For adsorbed particles on a surface, the model of an ideal adsorbed layer implies a lattice of 
places, and a multicomponent gas with components A = Z,Ai,...,A n where Z ia a free place 
and Ai are adsorbed particles on their places (each adsorbed particle occupies a place). There is a 
conservation law: 



i=0 



9 = const 



therefore, c = — Yll=i c * an( ^ me f ree energy has the form of the energy of the Fermi-gas: 



F(T,a,N) =F (T)a + Y,NMT) 

n 



1=1 



( ^) - 5i - 



T 



1 duijr) 
Rt dr 



dr 



(3.13) 



i=i 



i=i 



a 



where a is the surface area, F (T)a is the free energy of empty surface. 

For the systems, distributed in space, the density of free energy may be expressed through the 
concentrations, for example, for the perfect system (13.101) 



f(T,c) 



F(T,V,cV) 
V 

J2cMT) + RTJ2> 



In d 



T Q 



1 duj(r) 
Rt dr 



dr 



(3.14) 
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Therefore, for non-uniform system 

F= I f(T(x),c(x))dx. (3.15) 
Jv 

This formula is applicable if the space gradients are not too sharp. If f(T, c) is a convex func- 
tion of c then the minimizer of F (13.151 ) under given T, V and N = J v c(x dx (i.e. the equilibrium) 
is a uniform distribution and we should not expect spontaneous appearance of singularities on the 
way to the equilibrium. If f(T, c) is not a convex function of c then the minimizer of F may be 
nonuniform, non-smooth and non-unique: phase transitions are possible. In that case, the simple 
integral of the density / should be regularized by addition terms. Now, the standard approach gives 
the Ginzburg-Landau free energy: 



F= ip(T(x), c(x), Vc(x)) dx . 



v 



^(T(x), c(x), Vc(x)) = f(T(x), c(x)) + \ A,(V 



■i) 2 



(3.16) 



More general dependencies of Vc are also under consideration [58]. The chemical potentials fa 
for the free energy (13.161) are defined as variational derivatives of this functional, 

5F df(T,c) 

fa = t— = — t; A* Ac*. (3.17) 

Special analysis of the most general form of diffusion equations which provide the proper 
decrease of the Ginzburg-Landau free energy (13.161 ) was provided by Gurtin ||58i He introduced 
general nonlinear mobility matrices. 

The kinetic laws should satisfy the thermodynamic restrictions. That is, the dissipation should 
be positive. There are two physical forms of this law: (i) the available work should decrease and 
(ii) the entropy of the minimal isolated system, which includes the system under consideration, 
should increase. These two equivalent formulations correspond to two series of thermodynamic 
potentials: energetic or entropic series. There are several approaches to a general formalisms, 
which pretend to describe all systems that satisfy these monotonicity conditions ll56l [871 . Such 
approaches form the special discipline, nonequilibrium thermodynamics ll20l l59l l72ll or beyond 
equilibrium thermodynamics ||85ll . 

Our goal is different: we construct a method for assembling of a complex transport process 
from a mechanism combined by simple elementary processes. This approach for physics should 
satisfy the basic thermodynamic requirements. 



3.1.2. Markovian Microkinetics and Generalized Mass Action Law 

To satisfy the thermodynamic restrictions the kinetic law of the elementary reactions should have 
a special form, and the reaction rate constants for different elementary reactions should be harmo- 
nized. 
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Let us start from a reaction mechanism given by the stoichiometric equations (12.11 ): 



(3.18) 



where r is a reaction number, a ri and /3 ri are nonnegative numbers, the stoichiometric coefficients. 

In this Subsection, we return to some ideas of Michaelis and Menten [80] and to the Stueckel- 
berg analysis of the Boltzmann equation [97 ] and represent the general kinetic law for elementary 
processes. Detailed analysis was provided recently in |[54l . 

This law could be proved under the following assumptions: 

1. The elementary processes go through intermediate states (complexes or compounds) 11.441 



where p is the elementary process number; 

2. The amount of each compound B p is small enough to apply the perfect free energy formula 
(13.111) for them; 

3. The equilibrium between each compound and the corresponding linear combinations of 
reagents is fast enough to apply the quasiequilibrium approximation |[54l : 

4. The transitions between compounds could be described by a continuous time Markov chain 
(the Master equation or the monomolecular kinetics). 

The first three items of these assumptions correspond exactly to the celebrated Michaelis- 
Menten work ||80"1 . Later, Briggs and Haldane [7] abandoned the assumption of fast equilibria 
and produced the so-called Michaelis and Menten kinetics approximation. Original approach of 
Michaelis and Menten (fast equilibria with intermediate compounds + small amounts of com- 
pounds) was discovered again by Stueckelberg 11971 almost 30 years later. It was applied not to the 
kinetics of catalytic reactions but to the collision in the gas kinetics, as an alternative background of 
the Boltzmann equation. Gorban and Shahzad Il54l provided the detailed analysis of this approach 
to the general kinetic equation. 

Let the concentration of the intermediate compound B p (11.441) be <j p . The free energy (13.1 II) for 
the small admixture of compounds B p to the components Aj may be represented in the form 



(Fig.©: 




(3.19) 




(3.20) 



where q p is the standard equilibrium concentrations of B p . 

We assume that the standard equilibrium concentrations ?*(c, T) as well as the current concen 
trations q p are much smaller than the concentrations of Aj. 
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To formulate the results of Michaelis-Menten-Stueckelberg-Gorban-Shahzad (MMSGS kinet- 
ics) we have to introduce the basic notions in more detail. First of all, some of the formal linear 
combinations 

(a r , A) = a ri Ai , (/3 r ,A) = >^ (3 ri Ai 

i i 

may coincide. The same combination may be, simultaneously, the input combination of several 
reactions and the output combination of several other reactions. 

We assume that a fast intermediate compound B... corresponds not to a reaction but to a formal 
complex of the form (a r , A) or (/3 r , A) and this compound is the same for all reactions which 
include this complex. 

Let us call the formal linear combinations of the form (a r , A) or (/3 r , A) the complexes and 
enumerate them: l5 6 2 , ... , Q q . For each complex Qj, the corresponding vector of coefficients 
(a r or p r ) is yf Qj = A). 

The reaction mechanism (13.181) may be represented in the form Qj — > Q s for some pairs (J, s). 

The additional component, the fast compound Bj, corresponds to each complex Qj and the 
reaction mechanism Qj — > Q s (for some pairs (j, s)) (13.181 ) is extended to 

Qj ^ Bj — > B s ^ 6 S for some pairs (j, s) . (3.21) 

The fast equilibrium Qj ^ Bj gives 

i 

or 

9 = S' ( c ' T ) ex P ( ] > ( 3 - 23 ) 

where 

a/(c,r) 

is the chemical potential of A± and 

* = ln (f ) 

{RTdj = is the chemical potential of .&,). 

For the systems with fixed volume, the stoichiometric conservation laws of the monomolecular 
system of reaction between compounds are sums of the concentrations of Bj which belong to 
various connected components of the reaction graph. Under the hypothesis of weak reversibility 
there is no other conservation law. 

Let the graph of reactions Bj — > B% have d connected components C s and let V s be the set of 
indexes of those Bj which belong to C s : Bj G C s if and only if j G V s . For each C s there exists a 
stoichiometric conservation law 

& = J>- (3-24) 



60 



A.N. Gorban et al. 



Quasichemical Models of Multicomponent Diffusion 



For any set of positive values of /3 S (s = 1, . . . , d) and given c, T there exists a unique con- 
ditional maximizer <^ cq of the free energy (13.201 ): for the compound Bj from the sth connected 
component (j £ V s ) this equilibrium concentration is 

^ = R ^ (3 25) 

Inversely, positive values of concentrations ^ are the equilibrium concentrations (13.251) for 
some values of 8 S if and only if for any s = 1, . . . , d 

^■ = 0, if j,lEV s (3.26) 

{■dj = ln(<^/<^*)). This system of equations together with the equilibrium conditions (13.231) con- 
stitute the equilibrium of the systems. All the equilibria form a linear subspace in the space with 
coordinates fa/RT (i = 1, . . . , n) and dj (J — 1, . . . , q). 

Our expression for the free energy does not assume anything special about the free energy of 
the mixture of For the compounds Bi, we assume that this is a very small addition to the 
mixture of A i5 neglect all quadratic terms in concentrations of Bi and use the free energy of the 
perfect systems for this small admixture. 

The Master Equation for the concentration of Bj gives: 

-fo = Yl ( K i lQ ~ Kl ^) ■ ( 3 - 27 ) 

This kinetic equation should respect thermodynamics. For the Master equation this means the 
equilibrium condition 

Under this condition, the Master Equation (13.271) has the equivalent form: 

(!-!)■ 

In this form, it is obvious that <;* is equilibrium for the kinetic equation. 

All these expressions for concentrations of compounds and the Markov reaction rates result in 
the following kinetic law: 

1. The reaction rate of the reaction 0j — > Q s is 

w sj = (f) sj exp I J J , (3.30) 

where the quantity <p S j > is a kinetic factor, in the Markov model it corresponds to K S j^; 
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2. The kinetic factors (f) s j satisfy the complex balance condition: 



Yl fa* = Yl ^ for a11 • 



(3.31) 



in the Markov model this identity corresponds to the equilibrium condition (13 .28b ■ 
This is the macroscopic MMSGS kinetics. 

For this kinetics, the free energy decreases in time. To demonstrate this dissipation inequal- 
ity, let us formulate the macroscopic MMSGS kinetic equations in the original notations for the 
reaction mechanism (13.181) . Formula (13.301) for the reaction rate gives 



w r = (f) r exp 



RT 



(3.32) 



The kinetic equation under the isochoric conditions (the constant volume) are 

dc x - 

dl = ^ 7 ^' 



(3.33) 



where the stoichiometric vector 7 r = (3 r — a r . According to this equation, the dissipation rate is 

^ = V J>, lr )w r = VJ2 <fi r (», (A- ~ 0^)) exp , (3-34) 

where ^ = dFj dNi is the chemical potential. 

Let us introduce an auxiliary function like we did it for MAL (12.171) : 



0(A) = Y^ 4>r exp 

r 

This function is convenient because 

d0(A) 



\(a r ,/j) + (1 - A)(/3 r ,^) 
RT 



dX 



J^0 r (/i, (/3 r - a r ))exp 



A(a r ,/z) + (1 - A) 
RT 



Therefore, for the dissipation rate we get 



dF 
~dt 



-ve'd) 



(3.35) 



(3.36) 



(3.37) 



The function 9(X) is a sum of exponents. It is convex (9"(X) > 0). Therefore, if #(0) = 9(1) 
then 9'(1) > 0. This means that the identity 



9(0) = 9(1) 



(3.38) 
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is a sufficient condition for the dissipation inequality 

dF 

— < 0. 

dt ~ 

Some of vectors ot r , (3 r may coincide. Let there be q different vectors among them. We denote 
them y 1: ... , y q . For each y { we define Rf = {r\a r = yi}, R~ = {r \ f3 r = yi}. The sufficient 
condition for the identity (13.381) is 

(f> r =^2<f>r for all i . (3.39) 

reR+ reR- 

This condition is also necessary if we can vary <p r and fj^ independently and the Jacobian \dfXi/dcj\ 
has the full rank. This condition is the complex balance condition. 

Of course, the important particular case of the complex balance conditions is the detailed bal- 
ance condition: 

<t>t = <K , (3-40) 

where </>+, <p~ are the kinetic factors of the mutually reverse reactions: 0+ for (a r , A) — > ((3 r , A) 

and 0~ for (/3 r , A) — > (a r , A). 

3.2. From Cell- Jump Models to Continuous Diffusion Equations 
for Generalized Mass Action Law 

Let us start again from the stoichiometric mechanism of the diffusion: 

alA + £ o%A l l -+ Pli A l + E Pr\A? ■ (3.41) 

i i i i 

All the elementary processes have different numbers r and the space-inverted and reverse processes 
are represented separately. Let us consider all pairs of vectors, (a\, a") and (01, ffi) and numerate 
all the different pairs: yi, y 2 , ■ ■ ., where y q = (y l , y l q ). For each y q there are two sets of reactions, 

R + , R~: 

Rt = {r\ H, a?) = y q } , R~ = {r \ fi) = y q } . 
The reaction rate for the elementary process (13.411) is (13.321) : 

Wr = A exp (lMT)) + lo?M<?,T))j (3 42) 

The kinetic factors are functions of c 1 , c 11 and T. They should satisfy the identity of complex 
balance (13391) . 
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The cell-jump model gives us 

J = — •y r w r (c(x), c(x + /)) 



a*uv(c(x), c(x + /)) — ^ f3 l T w T (c(x), c(x + /)) 

SZ^U W r-(c(x),c(x + /)) - ^ W r (c(x),c(x + /)) 



The zero order term in / is 



It vanishes because of the complex balance condition. 
The first order term gives 



cxp 



RT 



)(E*(«?.v(£) 



Each term in this sum consists of the positive scalar pre-factor 



l exp 
and the matrix 



RT 



y\i \ ^a" ~~ ^"S 

multiplied by 



(3.43) 



(3.44) 



\RTJ 

This structure of the formula for the flux is very similar to (12.551) . Let us follow the same logic 
as in Subsection [23] to find the more convenient form of the expression for the flux (13.441) . First 
of all, let us group all terms with the same pre-factor. 
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For the set of all vectors y l q + y q l we use notation Z. For each z G Z 

Y z = {y q \y\ + y q l = z}. 

Analogously to (12.561) we get 



J = I / exp | 

zez \ 



RT 



■i 



Vq&z 3 \r&R+ reRa 



(3.45) 



Let us analyze J z for given z. 
First of all, 



Vq&fx 3 \r£R+ reRg / 







because each (p r from this sum enters it twice, with the opposite signs, one time as an element of 
the sum over R+ for some y q G Y z (with +) and the second time (with — ) as an element of the sum 
over R~ for another y s G Y z with the same sum 

yl + yf = y\ + yf = % ■ 

(Let us recall that 

a] + a 1 } = fi + Pi 1 

for all s and, hence, if r G R+ and al + a" = z then (3 l s + 0f = z and G 3^.) 

Let us mention that y 1 = z — yl 1 . Therefore, 



; '=-E«?E E - E ^ I v (w) • (3-46) 



In this formula, all the stoichiometric vectors are for the same cell, for the second one. 

Let us demonstrate similarity of the expression (13.461) to some formulas from the theory of 
Markov chains. Let us follow Subsection 12.5.1 and numerate the states of the auxiliary Markov 
chain by vectors y q G Y z . Elementary processes correspond to transitions between states. Each 
nonzero kinetic factor <p r corresponds to two vectors y q , y s G Y z : r G R q and r G R~ . 

Let us substitute the index r by two indexes q, s and use notation <f) sq for transitions y q — > y s 
(i.e. for the case r G Rt and r G R~). 

The complex balance condition (13.391 ) reads: 

$>«/ = E^- (3 - 47) 
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This is precisely the steady state condition for the Markov chain. 
Inequality (12.641 ) holds for the kinetic factors: 

E Kqi&qsKs ~ (psqKq) = E ^qsi^a ~ Tq) < (3.48) 

sg sq 

for any set of numbers n q . 

For the proof of the dissipation inequality it is convenient to rewrite J z in these notations: 

J Z =H Vl l E T.M ~ ■ (3-49) 

Vq&z j S 

Free energy in a domain V is 

F = / f(c,T)dx. 
Jv 

Let us estimate F in a bounded domain V with smooth boundary and without fluxes through 
boundary for isothermal conditions. 



dF 



Due to the representation of J z (13.491 ) 



qs 



(3.50) 



E ( v * (H) ' J /) = - ^ ' < 3 - 51) 



where n q is a space vector:^ = ^\ y%V x (j^) and for y q , y s 

vis + y l ls = z ■ 

This expression has exactly the form (13.481) for each space coordinate. Finally, 

EM-^U;)<o (3.52) 



\RTJ ' 3 

3 

and it is zero if all ix q = £\ V In q coincide. (The reverse statement is correct if the auxiliary 
Markov chain with the transition coefficients (p qs is ergodic for given z.) 
Therefore, 

F < 

for the generalized MAL (13.301) with the complex balance conditions (13.311) . 
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3.3. Detailed Balance for the Generalized Mass Action Law 
and the Dissipation Inequality 

3.3.1. Isothermal conditions 

Systems with detailed balance form an important subclass of the generalized MAL systems. 

Let us join each elementary process with its reverse process and represent the mechanism of 
diffusion by the pairs of mutually reverse processes: 

E <W + E ^ E &4 + E «a ? • (3-53) 

i i i i 

All the quantities for the direct process we mark by the upper index + and for the reverse process 
by the upper index _ . The simple algebraic relations hold (see also (12.4.1.1) : 

ajf * = l± and 7*f T = -7^ . (3.54) 

Therefore, 

<# + «" + = + fiT = <*l7 + = fit + • (3-55) 
Due to the principle of detailed balance, 



The reaction rates are 



,,1 „ „n 

I Pi . „ II 



wt(c\c 11 ) = r exp ( E a «^ + E a ™; 



,1 ,,n 



(c I , C II ) = r exp(^/3^ + E^ 



The cell-jump model gives us (13.431) 

J = E J r = ~ E 7r«(cF, C 11 ) - «, p -(cF, C 11 )) 



E 7^ [exp(E«H^ + E ) (3.57) 

\ i i 



exp 



3 I /fj 1 \ ^ oil /fj 

'RT ^ Pri RT 



In the first order approximation, we get analogously to (12.411) : 



J ri = -l(f) r exp ( E( Q; rj + a rj)^Fj E 7 «7rjV ( ^- 



RT / ^ \ RT 

3 / 3 



J r = — /0 r exp 



{a\ + aj. 1 ,//) 
RT 
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This expression for J r has the form of the multicomponent Teorell formula with the symmetric 
matrix of coefficients. This symmetry for nonlinear diffusion gives us the generalization of the 
Onsager reciprocal symmetry Il83ll84ll . We represented the nonlinear multicomponent diffusion a 
sum of elementary processes. For each elementary process 

where the scalar coefficient 

d r = Iffir exp 

is, from the thermodynamic point of view, an almost arbitrary positive quantity (because it includes 
the kinetic factor r ). "Almost" here means that some conditions of zero values (and the orders of 
these zeros) at the boundary when some of q = are prescribed by the factor 

and the logarithmic singularity of /!« when q — > 0. 

The internal symmetry of this formula makes the dissipation inequality obvious: in a bounded 
domain V with smooth boundary and without fluxes through boundary for isothermal conditions 

dF 




3.3.2. Generalization: Non-isothermal Processes 

Extension of the generalized MAL (13.301) on the non-isothermal processes is quite simple. Let us 
follow the paper [10J and include the "energetic component" Ajj in the list of components. Instead 
of the stoichiometric equations (12.11) . (13.181) we get: 

^2 a riAi + a Q Au ->■ ^2 PriA + /3qAu , (3.61) 
i i 

The corresponding macroscopic extensive variable for A v is the internal energy U with the density 
("concentration") u. To consider energy as the additional extensive variable, we should take the 
main thermodynamic potential for the isolated system from the entropic series. This is the entropy 

S: 

d S = -dU+-dV-Y^dN l . 

rp rp / j rp I 
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Let us extend the formulas for the generalized kinetics by additional component and take — -^ 
instead of j^. All the formulas including the dissipation inequalities remain the same. 
In isolated (isochoric) systems, U — and 



Ni 



IriWr 



where 

w r = <p r ex.p (a r , 

For transport processes, conservation of energy gives the following relations: 



a 



The space gradient of — enters the multicomponent Teorell formula as an additional force 
and the gradients of -j^ also enter the formulas for the heat flux. 
In particular, the simplest mechanism of transport, 



a Q A U ^ a Q A U 



generates the Fourier law: 



J Q = -la Q 4>exp (--^) V (- 



RT(x) 



-A(T)VT. 



The thermodynamic consideration cannot produce the temperature dependence of the thermal 
conductivity A(T) > 0. From the thermodynamic point of view, <p here is an arbitrary positive 
quantity. The problem of temperature dependence of A and its relations with other constants like 
diffusivity is widely discussed from the kinetic point of view [17711 . For computing thermal conduc- 
tivity various methods were developed including direct simulation and the Green-Kubo approach 
11631 . These methods were compared in 11931 . 

Thermodynamics may give the relations between different coefficients. For example, the prin- 
ciple of detailed balance produces the multicomponent Teorell formula with the symmetric matrix 
of coefficients (13.581) . The nonlinear reciprocal relations (13.591 ) could be automatically extended 
to the heat flux: just use the heat flux Jq as additional flux and —1/RT instead fXi/RT for the 
component A v . 

More relations we get for the specific mechanisms of transport. For example, for the activation 
mechanism of diffusion 



A 1 + a Q A\j ^ A u + a Q A% 



(3.62) 



the fluxes are: 



Ja 
Jo 



exp 
-/0exp 



RT 

A* - 
RT 



v(^) +ao v 



RT 



a «V(^)+aJV 



(3.63) 
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or in the matrix form 

For the mechanism (13.621) . the heat flux Jq is proportional to the diffusion flux J a with the 
coefficient oiq, that is the activation heat. In this activation mechanism, the activation heat travels 
with the particle from the cell I to the cell II. We can, for example, assume different behavior of 
the activation heat: let (3q distribute symmetrically after the diffusion jump: 

I i 

A 1 + a Q A\j ^ A 11 + -a Q A\j + -a Q A% . (3.65) 
For this mechanism, ^ v = — \ctu and 

(t)-^(^)d;t:)(vS,)' 

The heat flux Jq should be supplemented by the heat transport together with particles, £V /ij Jj, 
fTOl . The total heat flux is 

q = J V = Jq + ^ t^ J i ■ 0.61) 

i 

To describe the energy balance properly we have to include the work of various forces. The proper 
framework for modeling of the energy transport gives continuum mechanics. In its simplest form, 
fluid mechanics, we present these equation in the next Section. 



3.4. Momentum and Center of Mass Conservation 
3.4.1. Mass Transfer and Heat Transfer 

In this subsection, we briefly discuss coupling of the diffusion and thermal conductivity with fluid 
dynamics. The heat and mass transfer should satisfy the general laws of mechanics and, in partic- 
ular, does not violate the Newton laws. The diffusion and heat transfer equations do not present 
the complete theory and should be included into the context of continuum mechanics. 

First of all, let us introduce the mass average velocity. Let be the mass of mole (gram- 
molecule) for the component Aj. For each diffusion flux Jj the associated flux of mass is m 8 Jj. We 
introduced the fluxes Jj with respect to a frame, which is connected to our cells. For continuum 
motion, this frame should also move and we have to introduce the velocity of the frame, w. The 
flux of Ai associated with w is c,w. The corresponding flux of mass is mjQw. The total flux of Ai 
caused by diffusion and the frame motion is 

Ji = Ji + CjW . 

The mass density is 

P = miCi ' 
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the momentum density is 
the average mass velocity is 



y~]mjji 



v _ Tn^iJi =w | T^imJi 

Both the frame velocity w and the average mass velocity w are unknown. They are connected 
by the simple identity 

V = w + =^ , 

where the individual diffusion fluxes Jj are given by the mechanism of diffusion. 

The standard definition of the diffusion fluxes includes fluxes in the local center of mass frame 
where the average mass velocity is zero. Therefore, let us introduce the "proper fluxes": 

^ ^ ■ Tfl{ J % ^ ^ . TYl^ 

Cf i Ji VQ J i C{ ^ — , 3 % C{ _ — , . 

l^i m iCi m i°i 

These fluxes do not depend on the frame velocity. They are not independent and are connected by 
the momentum conservation: the total momentum is zero, 

^ rriiJi = . 

The heat flux in the local center of mass system is 



where Jq is given by the transport processes mechanism. 
For the energy flux, the standard approach 117011 gives 

fpv 2 \ 
v h u + P + Ju + viscosity terms . 



2 

The conservation laws give: 



d t p + div(pv) = , 

dt(pv) + div(pv <g> v) + VP = diva 



( ^- + u ) + div 



2 

d t Ci + div(vci + Ji) = . 



a : Vv 



(3.68) 



71 



A.N. Gorban et al. 



Quasichemical Models of Multicomponent Diffusion 



Here, a is the viscous stress tensor and o : Vv = £\. <?ij{dvi/dxj). The pressure P should be 
defined in accordance with thermodynamic properties of the mixture, for example, 

i 

where f(c, T) is the density of the free energy. 

The viscous stress tensor should be derived in (13.681) from the additional closure assumption. 
For the Newtonian liquid, 

^=H^ + ^J "v^J divv ' 

where fi (here and only here) is the shear viscosity and £ is the bulk viscosity. 

The individual equations in (13.681) are not independent. For example, the sum of the equations 
for conservation of Ni with coefficients gives us the first equation, the conservation of mass. 

The elastic energy and the various viscoelastic terms may be added to this picture. This is 
necessary to do and it is in our future plans. 



3.4.2. Mechanisms of Transport and the General Forms of Macroscopic Equation 

We developed a formalism of the mechanism of diffusion and heat conduction represented by the 
system of stoichiometric equations with the simple kinetic law exp(o:, fi/RT) and the balance 
condition (complex balance for the general Markov microscopic kinetics and detailed balance for 
systems with microreversibility). This formalism produces equations which are particular cases of 
the general nonequilibrium thermodynamic equations [|20l[85l . It is a very simple task to demon- 
strate that our transport equations are particular cases of the GENERIC formalism |[56l[85]|. Due to 
this two-generator formalism, evolution of any smooth function A of the state variables x is given 
by 

A A 

— = {A,E}+[AS] 

where E and S are the total energy and entropy, and {•, •} and [■, ■] are Poisson and dissipative 
brackets, respectively. 

The formulas for fluxes produced in this Section have the form of dissipative brackets: 

where M is a symmetric positively semidefinite operator, "the friction matrix". 

The general form of the "dissipative brackets with constraints" in application to multicompo- 
nent diffusion was produced very recently ll86ll . The flux of the ith component Jj in that formalism 
was presented by formula (54) ll86ll : 



^=-E a &[v(^)+a;v(4) 
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Our formulas belong to this type and give particular expressions for coefficients A. 

In the paper [86] a precise comparison of this formula with the classical expressions from EUl 
was presented and the equivalence of these general forms was proven. Now, we can just refer to 
these results. 

In addition to the general form, our approach gives the possibility to build the model from 
the elementary processes. This construction also satisfies the "constrains" (conservation laws) of 
diffusion because these conservation laws are implemented in the algebra of the stoichiometric 
coefficients (I2.22I ). 

4. Conclusion 

Chemical kinetics gave rise to the very seminal approach of the modeling of processes. This is, the 
stoichiometric algebra supplemented by the simple kinetic law. The results of this invention are 
now applied in may areas of science, from particle physics to sociology. In our work we extend 
the area of applications to nonlinear multicomponent diffusion. 

We demonstrated, how the mechanism based approach to multicomponent diffusion can be 
included within the general thermodynamic framework and proved the corresponding dissipation 
inequalities. To satisfy the thermodynamic restrictions, the kinetic law of an elementary process 
cannot have an arbitrary form. For the general kinetic law exp(a, jj/RT) (the generalized Mass 
Action Law), additional conditions on the set of kinetic factors <p are needed. There are two main 
sets of these conditions. The historically first of them, the condition of detailed balance, follows 
from the special property of the underlying microscopic dynamics, microreversibility. It was used 
by Boltzmann for his equation and then by many researchers like Wegscheider 111131 . Einstein Il25ll . 
and Onsager 11831 . The second and more general property we now call "semi-detailed balance" 
or "complex balance" was proposed by Boltzmann in his discussion with Lorentz. The theoretic 
principles underneath these conditions were discovered by Stueckelberg in 1952 [|97ll in application 
to the Boltzmann equation. It received the name "complex balance" ten years later in the works of 
Horn and Jackson [|62l and Feinberg [30]. Recently, it was demonstrated how this condition can be 
deduced from the quasiequilibrium and quasi steady state approximations from Markov kinetics 

El. 

Explicit formulas for nonlinear multicomponent diffusion combined from elementary pro- 
cesses can help in the practice of modeling (and already helped B51 ). The cell-jump formalism 
gives an intuitively clear representation of the elementary transport processes and, at the same time, 
produces kinetic finite elements, a tool for numerical simulation. 

There remain many questions for the future work: 

• It is necessary to extend the experience of modeling of real systems; 

• The analysis of diffusion in solids should be properly coupled with the mechanics of solids. 
The detailed quantitative theory of the Kirkendall effect may be a proper challenge here; 

• The mechanism representation should be extended to viscosity and viscoelasticity; 
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• The kinetic finite elements approach should be modified and extended to include continuum 
mechanics, perhaps, by proper joining with the lattice Boltzmann models approach; 

• Possible stoichiometric mechanisms of anomalous diffusion should be studied. 
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